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MIE:  A FORTRAN  PROGRAM  FOR  COMPUTING  POWER  DEPOSITION  IN 
SPHERICAL  DIELECTRICS  THROUGH  APPLICATION  OF  MIE  THEORY 


PURPOSE 

MIE  is  a FORTRAN  IV  program  written  for  the  IBM  360/65  system,  to 
compute  the  absorbed  power  density  at  internal  points,  the  a' erage 
absorbed  power  density,  and  the  total  absorbed  power  inside  a homogeneous 
spherical  dielectric  that  is  immersed  in  an  electromagnetic  field.  Re- 
sults are  obtained  through  applying  the  Mie  theory  and  using  double- 
precision arithmetic.  Here  double-precision  numbers  have  approximately 
16.8  decimal  digits  and  an  exponent  range  of  -78  to  +75. 

This  program  was  originally  produced  to  provide  special  test  results 
for  a computer  program  designed  to  solve  the  problem  of  ascertaining  the 
power  deposition  inside  arbitrarily  shaped,  finitely  conducting  biologi- 
cal bodies  exposed  to  electromagnetic  radiation. 

The  knowledge  to  be  gleaned  from  this  effort  is  directly  related  to 
the  research  effort  of  the  Radiation  Sciences  Division  at  the  School  of 
Aerospace  Medicine.  Briefly,  here  studies  are  being  currently  conducted 
(1)  to  determine  the  radiofrequency  radiation-induced  effects  in  biolog- 
ical specimens,  (2)  to  seek  out  possible  hazards  to  personnel  in  a radio- 
frequency environment,  (3)  to  accurately  measure  and  to  determine  the 
distribution  of  energy  in  the  whole  biological  body  or  just  in  a particu- 
lar organ,  (4)  to  find  ways  to  reduce  any  potentially  adverse  action 
between  RF  emitters  and  cardiac  pacemaker  and  prosthetics,  (5)  to  extrap- 
olate response  to  radiation  from  the  test  animal  to  man  in  a meaningful 
manner,  and  (6)  to  contribute  in  the  design  of  realistic  safety  standards 
with  a solid  basis. 

To  benefit  users  of  this  report,  program  MIE  is  described  in  suffi- 
cient depth  to  facilitate  implementation  on  any  modern  computer  and  job 
setups.  Detailed  coverage  is  provided  of  the  mathematical  theory  and 
formulas,  structure  and  sequence  of  control  parameter  and  data  cards, 
output  format,  and  function  subprograms  and  subroutines.  Included  in 
the  appendixes  are  an  extensive  discussion  of  the  Mie  solution,  sample 
problems  with  associated  computer  results,  and  a listing  of  the  program. 


MATHEMATICAL  DESCRIPTION 

We  consider  the  plane  electromagnetic  wave  that  is  irradiating  a 
homogeneous  spherical  dielectric  to  be  propagating  in  the  positive 
z-direction,  and  the  electric  field  E to  be  linearly  polarized  in  the 
x-direction  (cf.  Fig.  1).  A system  of  Cartesian  coordinates  with  origin 
at  the  center  of  the  sphere  is  used.  Also,  the  medium  which  surrounds 
the  sphere  is  taken  as  free  space  (c  vacuum).  Thus  our  embedding 
medium  is  a nonconductor,  and  both  the  surrounding  medium  and  the  sphere 
are  nonmagnetic. 
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Figure  1.  Directional  approach  of  the  incident  wave. 


The  Gaussian  system  of  units  (nonrat ionalized  c.g.s.  units)  is 
used.  In  this  system  the  macroscopic  form  of  Maxwell's  equations  may  be 
written  as 


7*D 

VxH 


VxE 


4irp, 
4 it 


1 SO 

J + - ~ 
c — c 8t 


1 

c 3t 


(1) 

(2) 

(3) 


V’H  = 0 » 


(4) 


with  D 
P 
H 
J 
E 
c 
t 


electric  displacement, 
free-charge  density, 
magnetic-field  intensity, 
current  density, 
electric-field  intensity, 
velocity  of  light, 
time . 


For  our  homogeneous,  isotropic,  permeable,  conducting  dielectric, 
and  linear  case,  the  constitutive  (material)  relations  are 
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D = e E , (5) 

B = H , (6) 

J = o E . (7) 

The  symbol  e represents  the  dielectric  constant,  and  o the  conduc- 
tivity. In  these  relations  and  in  equations  1-4,  the  magnetic  permeabil- 
ity, p,  was  set  equal  to  one. 

Since  we  are  considering  time-harmonic  fields.  Maxwell's  equations 
2 and  3 take  the  simpler  time-free  form 

VxH  = im2kli  , (8) 

VxE  = -ikH  , (9) 


where 


, 2n  w 

k = r = c • 

2 4ttq 

m = e - i 

0) 


with  k = propagation  constant  or  wave  number  in  free  space, 

X = wavelength  in  free  space, 
u = circular  frequency, 

m = complex  refractive  index,  a constant  for  a homogeneous 
medium. 

In  a homogeneous  medium,  relative  to  a Cartesian  coordinate  system, 
each  rectangular  component  of  vectors  and  H satisfies  the  scalar- 
wave  Helmholtz  equation 

V20  + m2k2tl>  = 0.  (12) 

At  the  same  time,  vectors  E and  H satisfy  the  vector  equation 

72A  + m2k2A  » 0.  (13) 

2 

Recall  that  the  formula  for  evaluating  V A is 


V A ■ V(7«A)  - 7x(7xA). 
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(Details  showing  that  vector  E satisfies  eq.  13  can  be  found  in  Appendix  A, 
pp.  19-21.) 

Now  two  independent  vector  solutions  of  equation  13  can  be  expressed 
as  vectors 


M = Vx(r  4>),  (15) 

—b  — — 

mkN  = VxM , , (16) 

-ill 'P 

where  the  scalar  function  il>  is  a solution  of  equation  12.  (Details 
pertinent  to  eqs.  15  and  16  can  be  found  in  Appendix  A,  pp.  27-30.) 

Since  r is  a constant  vector,  the  following  relationship  holds: 

mkM , = VxN  (17) 

— v — b 

Consequently,  if  we  choose  two  independent  solutions,  u ar>d  v,  of 
equation  12  and  construct  the  field  vectors  M , , N^,  N^ , we  will 

find  that  equations  8 and  9 are  satisfied  by  the  vectors 


E = M + iN  , 
— —v  — u 


(18) 


H = m(-M  + iN  ).  (19) 

— — u — V 

Here  i is  the  complex  unit.  (Details  pertinent  to  eqs.  18  and  19  can 
be  found  in  Appendix  A,  pp.  30-31.) 


Convenience  dictates  that  our  problem  be  embedded  in  a spherical 
system  of  coordinates  r(radial),  9(scattering  or  colatitude),  and 
^(azimuthal).  In  this  system,  the  components  of  vectors  and 

(i)  = u or  v)  can  be  expressed  as 


v = o 
r U’ 


N = 


mk  2 
3r 


(rib)  + mkr’b. 


(20) 


M9  rsinG  Deb  ’ 


N0  mkr  3r39  ’ 


(21) 


M*  - -7^* 


N.  = 


mkrsin<)>  3r3<J> 


(rib) 


(22) 


i 
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(Details  pertinent  to  eqs.  20-22  can  be  found  in  Appendix  A,  pp. 32-34.) 

♦ 

For  our  incident  plane-wave  radiation  expressed  in  the  form 


= a^_Eoexp[-i(kz-wt)  ] , 


H = a_^Eoexp  [-i  (kz-wt ) ] , (24) 

where  Eq  is  the  wave  amplitude,  and  polarization  vectors  a^  and  a^ 

are  unit  vectors  directed  along  the  x,y-axes,  van  de  Hulst  (5,  p.  122) 
chooses  the  independent  scalar  solutions,  u and  v,  of  equation  12 
for  the  induced  wave  as 


u = E exp  (io)t  )cos<f>  y me  (-i)  — ; — ttt-  P (cos0)j  (mkr),  (25) 

o n n(n+l)  n Jn 

n=l 


3exp(iu,t)sin$  I mdn(-i)n  (cose) jn (mkr), 

n=l 


The  expansion  coefficients  c and  d are  determined  by  using  the 

n ^ n 

appropriate  boundary  conditions;  pn(cos8),  the  associated  Legendre 

function  of  order  1;  and  j^(mkr),  the  spherical  Bessel  function  of 

the  first  kind  and  order  n.  (Details  pertinent  to  eqs.  25  and  26 
and  the  corresponding  forms  for  the  incident  and  scattered  waves  can 
be  found  in  Appendix  A,  pp.  34,  36,  37,  49-61.) 

The  boundary  conditions  on  the  tangential  components  of  vectors  Is 
and  II,  when  applied  to  the  components  in  equations  20-22,  yield 
a system  of  algebraic  equations  in  the  unknowns  a , b , c , and  d 
(van  de  Hulst's  notation  [5,  p.  123]):  n n n 


^(x)  “ an4n (x)  “ mcn^n ^ = 


K(x)  ~ Vn(x)  " CnK(y)  = °’ 


Vx)  " Vn(x)  ~ Vn(y)  = °* 


where 


^r\(x)  " bn?n(x)  “ m dn'J'n(y)  = °’ 


= ka, 


y = mka. 


i|<  (z)  - Zj  (z ) 


in(z)  . zh^<Z>  . 


hn  (z)  = jn(z)  - inn(z) , 


Hn  (z)  = Jn+^(Z)  - 1?WZ)- 


Here  x = MIE-size  parameter, 

a = radius  of  the  sphere, 

m = complex  refractive  index,  a constant  for  a homogeneous 
medium, 

X = wavelength  in  free  space, 

Jn+Xj(z)  and  ^n+i,(z)  = Vessel  functions  of  half-odd-integer 

order,  of  the  first  and  second  kind 
(Neumann  function)  respectively, 

(2) 

Hn+^(z)  = Hankel  function  of  half-odd-integer  order,  of  the 

second  kind . 

The  prime  on  $ ( z ) or  J (z)  denotes  first  derivative  with 
n n 

respect  to  the  function  argument. 

(Details  pertinent  to  eqs.  27-36  can  be  found  in  Appendix  A,  pp.  35-45.) 

Elementary  algebraic  manipulations  of  equations  27-30  result  in 

the  solution  of  c and  d : 
n n 


^(x)cn(x)  - * (x)t;(x) 

2 sr  

11  ~ m^n(y)C^(x) 


8 


(38) 


ip’Wz  (x)  - ip  (x)c'(x) 

, n n n n 

d 

n mij/  (y)  cn  (x)  - <|>n  (y)  4^  (x) 

(Details  pertinent  to  eqs.  37  and  38  can  be  found  in  Appendix  A,  pp.  42-49.) 

Utilizing  equations  18,  20,  21,  22,  25,  and  26,  we  can  derive  the 
r,  8,  and  4>  components  for  the  induced  electric-field  vector: 


2 2 

E^  = iEQexp(iuJt)cos<j>(m  krU  + 2mUR  + m krURR) , (39) 

where 

CO 

u = J1cn("i)n  pi(cos9)jn(mkr)-  (40> 

UR  = I c(-i)n  P^cosejj^mkr),  (41) 

n=l 


13  RR 


= J. 

n=l 


c (-i) 
n 


2n+l 
n (n+1 ) 


(cosS) j"  (mkr) . 
n ■ n 


(42) 


Ee 


EQexp  (iwt ) cos))  f mU 


+ i (~UR  + mURR)  ] 


(43) 


where 


U = 


n=l 


dn(-i) 


2n+l 
n (n+1 ) 


(cos0) jn (mkr) 


(44) 


UR 


n=l 


cn(-i) 


2n+l 

n(n+l) 


Tn(cos0) jn(mkr) , 


(45) 


URR 


n=l 


cn(-i) 


2n+l 
n (n+1) 


x (cos6) j ' (mkr) . 


(46) 


E = -E  exp(iojt)sin<t>[mU  + i(~UR  + mURR)]  , 

O Kl* 


(47) 
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where 


u=  ^ dn(-i)0  !lSiy  Tn(cOS0)j-(mkr)’ 


UR  = y cn(-i)a  2a+^--r  n (cos9)  j (mkr), 
L n(n+l)  n n 


URR  = y c (-1)  7 IT  (cosQ)  3 (mkr) . 

, n n(n-f-l)  n n 

n=l 


Here  a single  or  double  prime  denotes  a first  or  second  derivative  with 
respect  to  mkr  and 

tt  (cos0)  = — r~r  Pd(cos0),  (51) 

n smQ  n 

t (cos0)  = P3(cos0).  (52) 

n d9  n 

For  a sphere  of  radius  a (cm),  the  absorbed  power  density  at  an 
internal  point,  the  average  absorbed  power  density,  and  the  total  ab- 
sorbed power  are  computed  by  means  of  the  following  formulas  (with  o 
and  in  electrostatic  units  of  the  nonrat ionalized  c.g.s.  system): 

Pd(r,9,6)  = 0.05o  E-E*  watts/m3,  (53) 

-6  r fa  f2*  2 

P^ot.  = 10  V Sin9d*drd0  watts’  (54> 

1 0 1 O ' 0 

P = 10'P  / [ (4/3)ira3]  watts/m3.  (55) 

avg  tot 

Here,  * denotes  the  complex  conjugate.  In  program  MIE  a numerical 
integration  product  rule,  based  on  three  m-point  Gauss-Legend  re  quad- 
rature formulas,  is  used  to  evaluate  the  triple  integral. 

To  complete  our  discussion,  it  seems  appropriate  to  consider  the 
formulas  used  in  generating  the  values  of  certain  functions . The  formulas 


P3  . (cos9)  = 2n*d  cos9  P*  (cos6)  - -— — — ■ P d (cos0) , 

n+1  n n n n-1 


(57) 


sine  — P1(cos0)  = n cos6  P1 (cose)  - (nd-l)P1  ,(0036), 
at)  n n n-1 


together  with 


P^cosO)  = sin6, 


(58) 


P2(cos0)  = 3 sin6 cose 


(59) 


are  used  to  generate  function  and  derivative  values  of  the  associated 
Legendre  functions. 

Special  limit  values  are  also  obtained  by 


Lim 

6-+0 


P (cos6)  . , . 

n n(n+l) 

sine  2 


(60) 


Lim 

6-tt 


P (cos6)  . , .n+l  , , . 

n (-1)  n(n+l) 

sinO  2 


(61) 


The  forward  recurrence  relation 


n .,(z)  + n (z)  » ~+-1-  n (z) 
n+l  n-1  z n 


(62) 


is  used  together  with  the  relations 


nQ(z) 


cosz 


(63) 


^(z)  = 


cosz  sinz 


(64) 


to  generate  values  of  the  Neumann  spherical  Bessel  functions.  The 
generating  process  is  terminated  at  order  N when  the  following 
termination  criterion 


n (z)  > STOPR 

n 1 — 


11 


(65) 


r 


is  met.  Here  STOPR  is  a number,  say  1.0D25.  The  user's  needs  will 
determine  whether  or  not  STOPR  should  retain  its  presently  assigned 
value.  Our  own  demands  were  satisfactorily  met  for  both  real  ( 2tt r / f ) 
and  complex  (2nmr/f)  arguments  of  the  Neumann  functions  for  parameter 
ranges:  7 <_  |m|  <_  100,  5 <_  r <_  25  cm,  and  20  <_  f <_  1000  MHz. 

The  backward  recurrence  relation 


, „ 2n+l  . , s . , s 

Jn-!(Z)  = ~T~  JnU)  " Wz) 


(66) 


in  combination  with  an  appropriate  starting  value  is  used  to  generate 

values  of  the  spherical  Bessel  functions  of  the  first  kind,  ) (z). 

n 

This  technique  of  using  the  backward  relation  in  place  of  the  forward 
relation  helps  to  avoid  stability  problems. 

Values  of  the  derivatives  of  j (z)  and  nn(z)  are  obtained  by 
using  the  formulas: 

= 2^1  [nJn-l(z)"  <67> 


dz 


T Vz)  = 


z 


(68) 


h nn(z)  = -dn  [nVi(z)-(n+1)Vi(z)]' 


(69) 


PROGRAM  DESCRIPTION 


Program  MIE  consists  of  a driver  routine,  four  subroutine  subprograms, 
and  one  function  subprogram.  A list  of  the  driver  routine  and  subprograms, 
along  with  a brief  explanation  of  their  use,  follows. 

Driver  routine: 


MAIN  - used  to  input/output  data,  to  compute  spherical 
Bessel  function  values  for  a real  argument,  to 


compute  series  expansion  coefficients  c 


and 


and  to  direct  the  course  of  the  calculations. 


d , 

n 


12 


Subroutine  subprograms: 

EVEC  - generates  the  complex  radial,  colatitude,  and 
azimuthal  components  of  the  electric-field 
vector,  E,  and  the  scalar  product  E*E*. 

TERM  - generates  the  nth  power  of  -i,  where  i is 
the  complex  unit. 

BESS  - generates  an  array  of  values  of  each  of  the 
spherical  Bessel  functions  j (mx)  and 
n^Cmx)  for  complex  argument  mx.  Array 

dimension  = 100. 

PL  - generates  an  array  of  values  of  the  associated 

Legendre  function,  P^cosB),  and  an  array  for  its 
d ^1 

first  derivative,  — P^CcosQ),  for  n varying 
from  one  to  a maximum  N.  Array  dimension  = 100. 

Subroutines  BESS  and  PL,  plus  the  algorithm  (in-line  coded  in  MAIN) 

for  computing  the  spherical  Bessel  functions  i (x)  and  n (x)  for 

n n 

real  x,  are  modified  versions  of  subroutines  found  in  program  SUP[1] — 
one  of  two  allied  programs  developed  by  D.  S.  Drumheller  and  D.  E.  Setzer 
♦ (1)  as  a research  tool  in  the  study  of  problems  related  to  electromagnetic 

transmission  through  atmospheric  aerosols. 

Function  subprogram: 

GAUSS3  - a product  rule  for  the  numerical  evaluation  of  the 
triple  integral  for  the  total  absorbed  power  within 
the  sphere.  A basic  m-point  Gauss-Lsgendre  quadrature 
formula  is  used.  User's  option  is  available  for  the 
selection  of  m,  number  of  weighting  coefficients, 
and  associated  points  for  Gaussian  quadrature,  from 
the  set  of  integers  {2,3,4,5,6,8,10,12,14}  for  each 
of  the  three  Gauss  quadrature  rules  used. 

Blank  COMMON  is  used  by  driver  routine  (MAIN)  and  subroutine  EVEC. 

The  list  of  the  arrays,  variables,  and  constants  stored  in  this  area  is 


CN  - 


DN  - 


array  of  complex  expansion  coefficients  for 
coefficients  for  components  of  electric-field 
vector  E inside  the  sphere.  Array  dimension 
= 100. 
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A.IR 


AM 

AMK 


CEX 


DP 


P 


ALAMDA 

D1 

EO 

PHI 

PIE 

RKR 

STOPR 

SIXTH 


arrav  of  spherical  Bessel  functions  j^Cx), 

n = (complex  x)  values.  Array  dimension 

= 100. 

complex  index  of  refraction  (m) . 

complex  index  cl  refraction  (m)  times  wave 
number  (k) . 

complex  time-variation  factor  [exp(iwt)]. 

complex  index  of  refraction  (m)  times  Mie-size 
parameter  (x) . 

array  of  spherical  Bessel  functions  n (x)  , 

n = 1,...,N  (real  x)  values  or  first 

derivative  values  of  associated  Legendre  functions 

p^(9),  n = 1,...,N  (real  6).  Array  dimension  = 100. 

array  of  spherical  Bessel  functions  j (x), 

n = 1 , . . . ,N  (real  x)  values  or  values  of 

p^(0),  n = 1,...,N  (real  0).  Array  dimension  = 100. 

wavelength  of  incident  electric  wave. 

twice  the  radial  distance  to  an  interior  point 
of  the  sphere. 

intensitv  (field  strength)  of  incident  electric  field. 

azimuthal  angle. 

the  number  3. 14159265353973D0. 

the  reciprocal  of  Mie-size  parameter. 

a test  quantity  for  the  termination  of  the  recurrence 

relation  for  determining  the  spherical  Bessel  functions 

n (mx) . It  is  equal  in  magnitude  to  1.0D25. 
n 

the  value  of  sinh. 


THETA  - the  colatitude  angle  (9)  of  an  interior  point  of 
the  sphere. 

NC  - maximum  order  of  spherical  Bessel  functions  n (mx) 
(complex  mx)  less  two. 


A single  labeled  common  area,  GAUSS,  is  used  by  the  driver  routine, 
MAIN,  and  function  subprogram,  GAUSS3,  for  values  of 


f)  - diameter  of  the  sphere. 


Ml  - number  of  points  to  be  used  in  an  m-point  Gauss- 
Legendre  quadrature  rule  to  evaluate  innermost 
integral  of  triple  integral  for  total  power. 

M2  - number  of  points  to  be  used  in  m-point  Gauss- 
Legendre  quadrature  rule  to  evaluate  middle 
integral  of  triple  integral  for  total  power. 

In  subroutine  BESS,  if  variable  N,  the  maximum  order  of  spherical 
Bessel  function  r\n(mx)  (complex  mx) , tests  <_  2 , the  error  message 


PROCESS  CAN  NOT  PROCEED  SINCE  NM2  = 0 FOR  Z = 


is  printed  out  and  the  computer  run  is  terminated. 

In  function  subprogram,  GAUSS3,  an  invalid  value  of  one  of  the 
parameters  Ml,  M2,  or  M3  institutes  an  error-message  printout: 


ERROR  IN  INTEGRATION  CONTROLS.  Ml  = . . . 

M2  = . . . M3  = . . . 

A numerical  value  of  zero  for  the  triple  integral,  equation  54,  is 
returned  to  driver  routine,  MAIN. 

Input  to  program  M1E  is  by  keypunched  cards.  There  are  two  basic 
input  cards  with  structure  and  sequential  order  as  follows: 


Card  No.  1 (control 
Columns:  1-10 

11-20 

21-30 

31-40 

41-50 

51-60 

61-65 


parameters) 

FREQ.  Frequency  in  MHz.  (E10.3) 

EPS.  Relative  dielectric  constant.  (E10.3) 

SIGMA1.  Conductivity  in  mho/meter.  (E10.3) 

E0.  Intensity  (field  strength)  of  incident 
electric  field  in  volt/meter.  (E10.3) 

D.  Diameter  of  sphere  in  cm.  (E10.3) 

TIME.  Time  in  sec.  (E10.3) 

NOC.  Number  of  cases.  (15) 


66-70  IOPT  = 0:  Average  power  density  and  total 
power  not  computed;  = 1:  otherwise. 

71-73  Ml.  Number  of  points  for  Gauss-Legendre 
quadrature  rule  applied  to  innermost 
integral  of  triple  integral  for  total 
power  deposited  in  the  sphere. 
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74-76  M2.  Same  definition  as  Ml  except  appli- 

cation is  to  middle  integral.  (13) 


E 


77-79  M3.  Same  definition  as  Ml  except  appli- 

cation is  to  outermost  integral.  (13) 

Cards  Nos.  2-(N0C+l)  (coordinate  data) 

Columns:  1-10  R.  Radial  spherical  coordinate  of  inter- 

ior point  of  sphere  in  cm.  Range:  R > 0.  (E10.3) 

11-20  THETAD.  Colatitude  spherical  coordinate  of 
interior  point  of  sphere  in  deg.  Range: 

0 £ THETAD  <_  180.  (E10.3) 

21-30  PHID.  Azimuthal  spherical  coordinate  of 
interior  point  of  sphere  in  deg.  Range: 

0 £ PHID  <_  360.  (E10.3) 

The  last  card  of  a single  data  set  must  be  a termination  card  with 
the  symbols  /*  punched  in  columns  1 and  2.  Also  program  MIE  can  handle 
multiple  data  sets.  Each  data  set  [Cards  1-(N0C+1)3  is  stacked  one 
behind  the  other,  with  the  last  card  in  the  complete  data  deck  a termi- 
nation card. 

Under  option  control  IOPT=0,  program  output  printouts  consist  of 
the  title 

DEPOSITION  OF  POWER  INSIDE  OF  A HOMOGENEOUS  SPHERE  IMMERSED  IN  AN 
ELECTROMAGNETIC  FIELD 

followed  by  such  information  as 

FREQUENCY  = . . . MHZ 

WAVELENGTH  = ...  CM 

CONDUCTIVITY  = . . . MHO/M 

RELATIVE  DIELECTRIC  CONST  = . . . 

DIAMETER  = ...  CM 

TIME  = ...  SEC 

REFRACTIVE  INDEX  = . . . 

THE  HIGHEST  ORDER  NEUMANN  FUNCTION  COMPUTED  IS  OF 
ORDER  = . . . 

SIZE  PARAMETER  = . . . 


~ , V . 
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RATIO  OF  THE  FIRST  CALCULATION  OF  THE  ZERO  ORDER  SPHERICAL 
BESSEL  FUNCTION  OF  ARGUMENT  X TO  THE  VALUE  OF  SIN(X)/X  = 


RATIO  OF  THE  FIRST  CALCULATION  OF  THE  ZERO  ORDER  SPHERICAL 
BESSEL  FUNCTION  OF  ARGUMENT  Z TO  THE  VALUE  OF  SIN(Z)/Z  = 


INTERNAL  POINT:  RADIUS  = ...  CM  THETA  = ...  DEG  PHI  = . . . DEG 

ABSORBED  POWER  DENSITY  = . . . WATTS/M**3 

INTERNAL  POINT:  RADIUS  = ...  CM  THETA  = ...  DEG  PHI  = . . . DEG 

ABSORBED  POWER  DENSITY  = . . . WATTS/M**3 


For  IOPT  = I,  the  printouts  consist  of  the  information  set  forth 
above  plus  two  additional  statements: 

TOTAL  ABSORBED  POWER  = . . . WATTS 

AVERAGE  ABSORBED  POWER  DENSITY  = . . . WATTS/M**3 

Option  IOPT  = 0 with  NOC  = 0 produces  a printout  similar  in 
format  to  that  presented  on  page  16  but  without  the  data  on  the  internal 
points;  whereas,  option  IOPT  = 1 with  NOC  = 0 yields  a like  printout 
plus  total  absorbed  power  and  average  absorbed  power  density  results 
formatted  as  in  Che  above  paragraph. 

When  using  the  IBM  360/65  system,  the  combined  compilation  and  link 
editing  times  of  program  MIE  is  about  0.30  minute.  Execution  time  for 
specific  problems  will  be  found  in  Appendix  B. 
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APPENDIX  A 


MATHEMATICAL  ANALYSIS  OF  THE  MIE  SOLUTION 


Consider  four  vector-valued  functions  (E,  fl,  D,  and  B)  defined  on 
3 

the  Cartesian  product  (R  x R)  of  space  and  time.  Functions  E and 
H denote  the  electric-  and  magnetic-field  intensities;  functions  D 
and  _B  denote  electric  displacement  and  magnetic  induction  respectively. 

Vector  E>  is  also  called  electric  induction  by  O'Rahilly  (3,  p.  35). 

These  functions  are  related  by  Maxwell's  equations 

1 ^ 

-x-  + c 3t  = °*  (A-l) 


VxH 


c 3t 


4uaE 


(A-2) 


where  o is  the  conductivity  of  the  medium  at  (x,y,z)  at  time  t,  and 
c is  the  velocity  of  light. 

We  assume  the  constitutive  relations 


B = pH,  (A-3) 

D = e E , (A-4) 

3 

and  divide  R into  two  regions  in  each  of  which  u (magnetic  permeabi- 
lity), e (dielectric  constant),  and  o (conductivity)  are  constant 
functions. 

We  take  the  curl  of  both  sides  of  equation  A-l  and  substitute  into 
equation  A-2  and  have 


V(V-E)  - V E + (VxB)  = 0, 


(A-5) 


which  in  turn  implies  that 


2 

2 3D 

V(V-E)  - V E + — y 

c 3t 


. 3E 

4n  po  — 


3F  " °- 


(A-6) 
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From  equation  A-6  we  see  that  V*E  = p/e  implies 


V(p/e)  - A*  if  if  + ^£-|=0 

c 3t  c 


(A- 7) 


V_(V»E)  " V E + ^ — 2 + 

c 3t 


where  p is  the  charge  density. 


Let  us  assume  that  the  E field  is 


E = E exp (iwt) . 
— o 


Then  equation  A-2  implies  that 


VxH  = (—  E + — E )exp(iwt). 
c — o c — o 


4tt  no  — _ n 
2 3t  ’ 


(A-10) 


From  equation  A-10  we  see  that  if  e and  a are  constant,  then 


V*  (VxH)  = (IC0--4lT— )exp(iuit)V.E  = 0. 
— c -o 


V • E =0 


(A-l 1) 


(A-12) 


in  a region  where  e and  o are  constant.  Therefore,  for  time- 
harmonic  waves,  equation  A-8  may  be  replaced  by 


? d L , do 

-7  E +££_!+  _Z 

--  c2  2 2 3t 

c 3t  c 


in  a region  where  e and  o are  constant.  Substituting 
equation  A-9  into  equation  A-13  we  find  that 


( A- 1 3 ) 
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i 


(A-14) 


V E + ( 


peu 


^)E  = 0. 


Let 


A2  = w2[^  - i = (— )2(pe-i  ^H£). 

2 2 C U) 

c wc 


(A-15) 


Then  substituting  equation  A-15  into  equation  A-14,  we  have,  with 
2 

k=u>/ c and  m = ye-i4tpa/w,  the  vector  Helmholtz  equation 


2 2 2 
V_  E + m k E = 0, 


(A-16) 


where  k and  m are  the  wave  number  and  complex  index  of  refraction, 
respectively,  and  remain  as  such  in  the  rest  of  the  paper. 

Now  the  scalar-wave  equation 


2 2 2 
Vi|/  + mkt|j=0 


(A-17) 


in  spherical  coordinates,  a separable  coordinate  system  for  the 
equation,  has  the  form 


3 , 23<k  , 1 3 , . . 3ik  . 1 • 3 * . 2,2, 

^(r  37>  + "TT7  30'(sine3e)  + T ~~2  “A  m * = °‘ 

n r sin  0 dp 


r sin0 


Let 


(A-18) 


<P  = R(r)0(0)*(<|>). 


(A-19) 


Substituting  equation  A-19  into  equation  A-18  yields 


04,^T  l?^2  R')l  + -2^—  fe-<sin0  ©' ) 

r r sin0 


+ ~2  4>"  + m2k2R  0 4>  = 0. 

r sin  0 


(A-20) 
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[ 


Expansion  of  terms  in  equation  A-20  results  in 


e*[R"  + (2/r)R']  + O’  + ©") 

r 


+ -„9  *v  4>"  + m2k2R04>  = 0, 


2 2 
r sin  0 


(A-21 ) 


and  dividing  all  terms  in  equation  A-21  by  R0<t>,  we  infer  that 


i[R"  + (2/ r )R’  ] +[^1(0"+  0')  + -T1- 


; L 2 sin0  ” ' ' 2.2.. 

r sin  04> 


2,  2 


+ m k = 0 . 


(A-22) 


Finally,  multiplying  all  terms  of  equation  A-22  by  r“~,  we  obtain 
the  separated  form  of  the  equation 


{^-[R"  + (2/r )R ' ] + m2k2r2}  + -|(0"  + gff  0') 


2 

sin  0$ 


<J>"  = 0. 


(A-23) 


We  use  the  following  lemma: 

herma  A-l.  Let  y be  a Bessel  function  of  order  (n  + *s) , then 


x2y"  + x y'  = [n(n+l)  + h - x2]y 


(A-24) 


along  with 


v * n(n+l) 


(A-25) 


and 


u =y//x~ 


(A-26) 
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* r 


implies  that 


x2u"  + 2xu'=  (v2- x2)u  . 


(A-27) 


Proof  of  Lemma  A-l.  Equation  A-24  is  satisfied  because 
2 2 

(n+^j)  = n"  + n + \ = n(n+l)  + \ and  by  the  definition  of 

Bessel's  equation  (Whittaker  [6,  p.  38]).  Differentiating  u 
we  find  that 


u'  = y ' / V"x  - (1  / 2 ) y / x3/2 , (A-28) 

u"  = y'vVjT"-  (l/2)y'/x3/2  - (l/2)y'/x3/2 

+ (3/4)y/x5/2,  (A-29) 

or 

u"  = y"/Vx”  - y'/x3/2  + (3/4)y/x572  . (A-30) 

Now 

x 2 u " + 2xu'  = x372y"  -\Cy'  + (3/4)y/V^~  + 2Vxy' 

- y/VT~  (A-31) 

or 

x2u"  + 2 xu'  = a/^)[x2y"  + xy'  - (l/4)y].  (A-32) 

From  equations  A-24  and  A-32  it  follows  that 


x2u"  + 2xu'  = (1/Vx7{[n(n+1)  + (1/4)  ]y  - (l/4)y-x2y}J 
or 


23 


x~u"  + 2xu'=[n(n+x)  _ x 2 ] (y  f\j x ) , 


=[n(n+l)  - x )u. 


(A-33) 


This  completes  the  proof  of  Lemma  A-l. 

2 

Now  equation  A-23  implies  that  there  is  a constant  v such  that 


(l/R)(r2R"  + 2rR')  + m2k2r2  = \ ? 


(A-34 ) 


r2R"  + 2rP ' + (m2k2r2-v2)R  = 0. 


(A-35) 


Let  us  write 


R(r ) = z (mkr) . 


(A-36) 


On  substituting  equation  A-36  into  equation  A-35,  we  obtain 


2 2 2 ? 7 9 2 

mkr  z" (mkr)  + 2mkrz'(mkr)  + (m  k r -v  ) z (mkr)  = 0, 


(A-37) 


which  according  to  Lemma  A-l  is  true  provided  that 


z(mkr)  = C 


Zn+^(mkr) 
~\J  mkr 


(A-38) 


where  C is  an  arbitrary  constant 

2 

v = n(n+l), 


(A-39) 


and  is  a Bessel  function  of  order  n-f-^. 

Substituting  equations  A-36  and  A-38  into  equation  A-23  results 
in  the  relation 


"<"«>  +i(0" + nff  • °- 

sin  6$ 


(A-40) 


Now  multiplying  all  terms  of  equation  A-40  by  sin^e,  we  find  that 
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[n(n+l)sin20  + (1/0)  (sin290M  + sin0 cos90' ) ] + (l/4>)4>"  = 0.  (A-41) 


(1/*)*"  = -t. 


(A-42 ) 


#"  + rt  = o 


(A-43) 


implies 


<*> ( <J> ) = c^sin(£(j>)  + c^cos  (•£$) 


(A-44) 


for  some  constants,  and  C2*  Substituting  equation  A-42  into 

equation  A-41  yields 


n(n+l)sin2O0  + sin290"  + sin0  cos00'  - £20  = 0 


(A-45) 


(l-cos20)0"  + sin6  cos90'  + [n(n+l)  (l-cos*^)-/^2  ]0  = 0.  (A-46) 


Let  us  set 


0 = w(cos9) , 


(A-47) 


0'  = -sin0  w'(cos0), 


(A-48) 


0"  = -cos0  w' (cos9)  + sin20  w"(cos9).  (A-49) 


Substituting  equations  A-48  and  A-49  into  equation  A-46,  we  obtain 


L 


(l-cos20) [-cos0  w'(cos0)  + sin"9  w"(cos0)] 


2 2 „2 

+ cos0  sin  0w'(cos0)  + [n(n+l) (1-cos  Q)-t  ]w(cos0)  = 0 


(A-50) 


and  dividing  all  terms  of  equation  A-50  by  (l-cos“6)  yields 


(1-cos  9)w"(cos0)  - 2 cos9  w' (cos0) 


+ [n(n+l)  - ]w(cos6)=  0. 

1-cos  0 


(A-51 ) 


Now  set  x = cos9  in  equation  A-51  and  discover  that  w(x) 
satisfies 


(l-x2)w" (x)-2xw' (x)  + [n(n+l)  — v]w(x)  = 0.  (A-52) 

1-x 


It  can  be  shown  that 


w(x)  = P^(x) 


(A-53) 


l 

is  a solution  of  equation  A-52,  where  Pn(x)  is  a solution  of  the 
associated  Legendre  ordinary  differential  equation  (Whittaker  [6,  p.  324]) 


Hence 


Zn+U(mkr)  p 

— P^  (cos0)  [c^cos  (£<S)  + C2sin (£<J>)  ] , 


V 


(A-54) 


mK.r 


where  denotes  a Bessel  function  of  order  n-f^.  Observe  that 


= R(r)0(0)4>(d>) 


(A-55) 


satisfies  equation  A-17. 
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1 


The  following  vector  functions: 


l = n. 

(A-56) 

M = Vx  (a_tp ) = Vpxa  , 

(A-57) 

N = (l/mk)V>;M, 

(A-58) 

in  which  a is  a constant  vector  and  p is  a spatial  function,  satisfy 
the  same  vector-wave  equation  since  the  components  of  these  vectors  are 
linear  combinations  of  derivatives  of  solutions  of  the  basic  scalar 
Helmnoltz  wave  equation — a linear  partial  differential  equation  with  con- 
stant coefficients. 

The  validity  of  equation  A-57  is  based  on  the  following  proposition: 

Proposition  A-l.  For  every  C1"  function  p and  every  constant 
vector  a_,  we  have 


Vipxa  = 


i 

i 

k 

_3£ 

dip 

8x 

3y 

dz 

ai 

a2 

a3 

. , 94>  . , dip  dip. 

= — a35V  " a2^)-i(a3^  ' al^ 


+ w,  M _ a li} 

+ — ^a23x  al3y) 


(A-61) 


Since  the  right  sides  of  equations  A-60  and  A-61  are  identical,  the 
proposition  is  proven. 

The  approach  of  Stratton  (4,  pp.  392-423)  is  to  express  the  vector 

potential  A as  a linear  combination  of  vector  functions  L , M 

— — n — n and 

N where  these  depend  on  ip  , the  nth  spherical  harmonic  associated 

with  the  scalar  Helmholtz  wave  equation. 

Observe  that 


V_x(ri|j)  = Vipxic  + ipVxr, 


(A-62) 


where  _r  is  the  position  vector.  To  simplify  equation  A-62,  we  need 
to  represent  the  basic  operators  in  an  orthogonal  coordinate  system. 


Let  (x,y,z)  denote  a point  in  Cartesian  coordinates.  Suppose 

that 


12  3 12  3 

x = x (u  , u , u ) , y = y(u  ,u  ,u  ), 


12  3 

z = z (u  ,u  ,u  ) , 


(A-63) 


If  the  position  vector,  _r,  is  expressed  as 


r = xi 


+ yj_  + zk  (i,i,k  unit  base  vectors) , 


(A-64) 


then  in  Cartesian  coordinates 
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Vi^xr  = 


ip  ip 

x y 


- i(zii  - y\p  ) - j (zif/  — xip  ) k (y^>  - xip  ) . 

y z x z — x y 


(A-65) 


It  is  therefore  natural  to  define  operators  L , L , and  L by  the 
rules  x V z 


V - z “ y If  ’ 


L * = x ! 1 _ z |i 

y 3z  3x 


L p = y - x ^ 
y 3x  3y 


(A-66) 

(A-67) 

(A-68) 


Note  that 

2 


l.»%  - *<%-  + 4 + ^V-)  - + A,  + 

3x  3y  3y  3z  3y  3z3x  3z3y  3z 


[(— )2  + (— )2  + (— )2](z  - x 

U3x'  ^3y;  (‘3z;  U 3y  3y;  * 


(A-69) 


2 2 

That  is  to  say,  L 7 = V L . Moreover,  we  have  the  following  theorem: 


Theorem  A-l.  For  all  C functions  ip,  we  have 


2 2 


7 L p = 

X 

L 7 ip  , 

X 

(A-70) 

2 

2 

7 L ip  = 

y 

L 7 ip  , 

(A-71) 

2 

2 

7 L ip  = 
z 

L 7 ip  . 
z 

(A-72) 

29 


Hence , if 


2 2 2 
V <P  - -m  k ip, 


the  vector  function 


satisfies 


V = iL^  + + kL^ 

= VlpXJT 


2 2 2 

V V = -m  k V. 

Proof  of  Theorem  A-l.  If  V is  defined  by  equation  A-74, 
then  equations  A-70  through  A-73  imply  that 


2 2 2 2 

7V  = iVLiJ)+iVLii  + kVLi|) 
- — x y — z 


2 2 2 

= iLVil/+iLVij;  + kLViJ; 
— x ^ y — z 


= + 2Ly  (-m2k2ii)  + kL^  (-m"k2il<) 


2,  2„ 

= -m  k V.  (A- 


This  establishes  Theorem  A-l. 

If  u and  v are  two  solutions  of  the  scalar-wave  equation 
Maxwell's  equations  are  satisfied  by 


H = My  + iN^  (i  is  the  complex  unit).  (A 

Equation  A-l  Implies 


VxE  + ^ b = 0. 


(A- 


(A-73) 


(A-74) 


(A-75) 


76) 


then 


77) 


78) 
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Thus 


H = 


i ( — - — ) VxE  . 

1IM 


Furthermore,  in  view  of  equations  A-58  and  A-77 


H = (— )mkN  + (— )[Vx(VxM  )/mk]i 
yw  — v yw  — u 


or,  equivalently. 


H = lmkc  N + 

- yw  -v  Vwmk'1-'-^'  -^uJ 


•)  [V(V«M  )-V  M ] i 


The  definition  of 


M 
— u 


by  the  rule 


M = Vxru, 
u 

the  vector  identity  (the  divergence  of  a curl  is  zero),  and 
A-75  or  A-76,  imply  that 


H = ^ N 

— UU)  — \i 


+ ^ M i 
yw  — u 


or 


H = 


mkc 

yoj 


(-M  + iN  ). 
— u -v 


In  the  Gaussian  units  for  a nonmagnetic  body,  y = 1. 
k = w/c,  we  have 


H = m(-M  + 
— — u 


iN  ) 
—v 


(A-79) 

(A-80) 

(A-81) 

(A-82) 

equations 

(A-83) 

(A-8h) 
Thus,  since 

(A-85) 
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In  the  spherical  coordinate  system  with  unit  base  vectors,  e , e 


an^  vector  A can  be  represented 


r’  -9, 


A = A e + A.eQ  + A e . 
r r 0 — 0 (j)—(p 


(A-85.1) 


Then 


-x-  = {[^(rsln9V-  ^rV^r 


r sin9 


3Ar  3 

+ lTf  ~ 87  (rsinQA^r^ 


3 3A 

+ [TF(rAe>  ' ¥lrsin6^: 


(A-86 ) 


We  take  Af  = r*,  A = 0,  and  A.  = 0,  and 


observe  that 


M = Me  + Me  + M e 
rr  9^0  cp — <|> 


(A-87) 


implies  that 


M = 0, 
r ’ 


M0  r sine  3(^r^  * 


(A-88) 

(A-89) 


M<!)  " r 90  • 


(A-90) 


Furthermore,  the  functions  N satisfy  the  relation 


\ - 5k  1% 


(A— 9 1 ) 
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' - ~ . v-  _ 


fc. 


which  implies  that  if 


N,  = N e + N e + e , , 

-tJ/  r~r  e-  9 <p4> 


(A-92) 


then 


Nr  = ^ {“  (7~'  - f 3?  [sin6  3?  <r*>  ]"  (r21'.  27}“2  (r*>  }* 

r sm6  r sin  9 3$ 

(A-93) 


Ne  = ^ {'  17  [rsin9  (-  7^(r^))]r}(-T-T) 

r sinB 


— [—  — — (rijO  ] 
mk  r 3r36  * J ’ 


(A-94) 


= \ Y1.  ■ ^(^)Hrsine. 

mkr  sin9 


(A-95) 


Hence 


N = = — — 

4>  mkr  sin9  3r3<{i 


(A-96) 


Now  equation  A-93  may  be  simplified  by  using  the  fact  that  tl 
satisfies  equation  A-18.  We  have 


N = i_  [II  (r2  li' 
r mk  lr  3r  3r 


) + 


2,  2 , , 

m k r ] 


( A-97 ) 


Since 


+ i, 


(A- 98) 
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partial  differentiation  with  respect  to  r yields 


2 2 

3 / o dJk  x 1 3*  1 d , 23J\ 

— (r*)  = r 72  + 2 3?  = 7 37(r  3?  ■ 


(A-99) 


3r 


3r 


Substituting  equation  A-99  into  equation  A-97,  we  see  that 

2 

Nr  = i [72(r^  + n’2k2r'1'1  • 

9r 


(A-100) 


The  incident  wave  is  described  by  the  vectors 


E = £xEQexp[-i(kz-(ot)  ] , (A-101) 

H = a E exp[-i(kz-ujt ) ] , (A-102) 

— — y o 


where  vectors  a and  a are  unit  vectors  directed  along  the 
—x  — y 

x,y-axes.  Vectors  _E  and  H can  be  expressed  also  in  the  form  of 
equations  A-77  and  A-85  with  scalar  functions 


u - E exp(iojt)cos<}>  7 (-1)  . . P (cos9)i  (kr)  , (A-103) 

o u , n (n+i ) n n 

n=l 


00  >■)  - i *i 

v = E exp(iwt)sin<|!  £ (-i)  ~ (cos9)  (kr ) 


n=l 


n(n+l)  n 


(A-104) 


The  representation  of  vector  E is  given  by 


E = M + iN  . 
— ^v  — u 


(A-105) 


The  terms 


” IKe  S - V9-*)v’ 


(A- 106) 
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2 


(A-107) 


(N  ) 


— u 


e 


_L  I 3 

mk  r 3r39 


(ru)  = Gg  (9  , (J>) 


make  up  E_. 

9 


Similarly,  for  E, 


(M  ), 
— v $ 


3v 
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F (6,<j>)v, 


o;  ).  = 


1 3_ 

mkr  sin8  3r3b 


(ru)  = 


( A— 1 08 ) 


(A- 109) 


This  does  not  remain  true  when  in  equations  A-103  and  A-104  the  (-i)n 

is  replaced  by  (-i)na  and  (-i)nb  , respectively,  where  the 

n n 

(a, , a_  , . . . ,a  , . . . ) and  (b. ,b„  , . . . ,b  , . . . ) are  members  of  the 
* 2 n 1 i.  n 

sequence  space  appropriate  to  the  functions  being  represented. 

We  see  by  differentiation  in  the  tangential  directions  that  if 

1 3 

- -(ru)  and  v are  continuous  in  r for  each  6 and  b,  then 

Ea  and  F.  are  continuous  across  the  boundary  of  the  sphere. 

9 0 

Now  we  repeat  this  argument  for  the  H vector.  The  H vector  is 
defined  by 


H = m(-M  + iN  ) . (A-110) 

To  get  the  tangential  components  of  H,  simply  interchange  u and 
v in  the  preceding  argument.  Thus,  the  continuity  of  tangential 

H is  assured  if  ^ and  mu  are cont inuous  in  r for  each  9 

and  b . 

We  introduce  a new  set  of  functions  which  differ  from  spheric  n 
Bessel  functions  by  an  additional  factor,  z.  These  functions  arc 
defined  by 

C>n(z)  = z j n (z ) = (ttz/2)  , ( A- 111) 

X (z)  = -zn  (z)  = " (nz/2)^N  , (z),  i A— 112) 

n n n+'i 
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(A-113) 


Cn(z)  = zhn(2)(z)  = (irz/2)V^(z), 


where 


Jv(z)  - 


-l)k(z/2) 2k 

k=0  r(k+i)r(k+v+i) 


[ l 


Hz/2)v 


(A— 1 14) 


is  holomorphic  for  integral  v and  converges  for  all  v if  the 
appropriate  definition  of  (z/2)v  is  given  and 


Y 

v 


(z) 


Jv(z)cos(vn)-J_^(z) 
sin  (vtt) 


(A-115) 


H<2)(z)  - Jv(2)  - 1 Yv(z) 

J_v(z)  - exp  (ivir)  (z) 

= -i  sin  (vr)  ’ (A-116) 

\(z)  = Yv(Z).  (A-117) 

(The  functions  are  called  Bessel  functions  of  the  second  kind, 

Neumann  functions, or  Weber  functions.) 

An  outgoing  scattered  wave  can  be  obtained  by  applying  a linear 
operator  to  the  ordered  pair  (u,v),  where 

00  o 

u = Eoexp(iwt)cos4>  ^-an(-i)n  (cos9) hn(2)  (kr ) , (A-118) 

00 

v = Eoexp(i(i,t)sin^  [ -^(-i)"  (cose)htJ2)  (kr) . (A-119) 

n»l 
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Similarly,  an  induced  wave  can  be  obtained  by  applying  a linear  operator 
to  (u,v) , where 


u = Eoexp(iwt)cos$  l mcn(-i)n  |Y^yPn(cos6)jn(mkr)  * (A-120) 

n=l 


v = Eoexp(iwt)sin<{>  £ md^C-i)11  n^+1)P^(cos6)in(mkr)  • (A-121) 


n=l 


The  boundary  conditions  are  defined  by 


incident  , -scattered  -inside 
E + E = E 


tangent 


tangent 


tangent 


(A-122) 


and 


incident  .scattered  inside 
H + H — H 


tangent 


tangent 


tangent 


(A-123) 


Since  the  tangent  space  of  the  sphere  is  two-dimensional,  we  have  four 
equations  in  four  unknowns.  These  are  assured  by  the  continuity  of 


1 3 (ru)  , 3 (rv) 

mu’  m v’  and 


Thus,  equations  A-103  and  A-104, 


equations  A-118  through  A-123,  and  the  continuity  relations  imply  that 
the  relations  that  must  be  satisfied  are 


(mu)  + (mu)  = (mu) 

(A-103)  (A-118)  (A-120) 


( A— 1 24 ) 


<1  JLaTi>  + «;  ifri)  = ^ • a-125) 

(A-103)  3 (A-118)  m 3r  (A—  3 ’ 


V(A-104)  + V (A-119)  V(A-121) ’ 


(A-126) 


+ m /ilrvix 

3r  (A-104)  3r  (A-119)  3r  (A-121) 


(.A-127  ) 
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In  equations  A-124  and  A-125,  the  in  denotes  the  refractive  index  of 
tb-  propagating  medium.  Outside  the  sphere,  m = m^  = 1 . 

Equation  A-124  with  r = a yields 

CO 

m E exp(iwt)cosij)  7 (-i)n  ~~ - P1(cos6)j  (ka) 
o o u n^n+1)  n n 

n=l 

CO 

+ m E exp (iojt)cos<{)  7 -a  (— i ) n “-7— 77.  P'*'(cos0)h^^  (ka) 

00  _i  n n(n+l)  n n 


mE  exp(iu)t)cos<(>  7 me  (-1)  — , , , P (cos9)j  (mka) , 

o n n(n+l)  n n 

n=i 


Multiplying  all  terms  of  equation  A-128  by 


setting 


x = ka. 


y = mka  = mx. 


and  using  the  fact  that 


m = 1 , 
o 

we  finally  arrive  at 

(2)  9 

xi^Cx)  - anxh^  (x)  = m“cnxjn(mka) 

= mcnyjn(y) 

= me  ijj  (v) . 
n n ‘ 

Equations  A-lll  through  A-113  and  A-132  yield 


(A-128) 


(A-129) 


(A— 1 30) 


(A-131) 


(A-13.) 


il»  (x)  - a c (x)  = me  i|>  (y). 
n n n n n 


(A— 1 33 ) 
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■V  r 


Next  we  make  use  of  the  continuity  of  — ^ru-^  . Using  equations 

m dr 

A-103,  A-118,  A-120,  and  A-125  results  in 


A vxp(iwt)cos^  u-v*  ^Si)pi(cose)  f?(krjn(kr)) 


+ Ak  Enex?  UwOcos*  l (-a  (-i)n)  — i-.P^(cose)  |-(krh(2)(kr)) 


IT.  K O 
O 


n(n+l)  nvw  ~ ' 3rv  n 


mk  EoeXp(i't)cos?  n=^mCn (-i)n  ^Sr)PAcos0>  (mkr ) ) . (A-134) 


Thus,  using  equations  A-lll  through  A-113  and  the  orthogonality  of 

the  P (cosS),  we  arrive  at 
n 


*n(x)  " Vn(x)  = cn*n(y)  ' 


(A-135) 


Here  we  have  used  the  fact  that  y = mkr  and  that 


— [mkrjn(mkr)]  = mk  [ yJn (y) ] 


(A-136) 


The  continuity  of  v is  expressed  in  equation  A-126,  and  use  of  A-104, 
A-119,  and  A-121  shows  that 


CO  Oil 

Eoexp(iut)sin$  ][  (-i)"  —^P^cosO) Jn<kr) 


+ E exp  (loot)  sin<i>  £ -b  (-i)n  2/+l1',P1(cose)hU)  (kr) 
o n n (n+1)  n n 


03 

EoexP(iwt)sin0  [^d^C-i)0  (cos0)  j^(mkr)  , (A-137) 


Thus,  multiplying  both  sides  of  equation  A-137  by  kr  and  using 
equations  A-lll  through  A-113,  A-129,  and  A-130,  we  have 


' 


'l'(x)  -bC  (x)  = d iJj  (y) . 
n n n n n 


(A-l 38) 


and 


T,  . r 3 (rv) 

The  continuity  of  — r — - 
3r 

this,  in  view  of  equations 


is  expressed  by  the  equation  A-l 27 
A-104,  A-119,  and  A-121,  implies  that 


E exp(iwt)sinb  j (~i)n  P1  (cos9)~[ rj  (kr)J 

o n (n+1 ) n dr  n 


+ E exp(iwt)sin<}>  [ -b  (-i)n  P1  (cos0)|-[ rh(2  5 (kr ) ] 

° n n(n+l)  n dr  n 


= Eoexp(iu,t)sin<>  l mdn(-i)n  (C0S9)|_[r  (mkr ) ] . 


n=  1 


(A- 139) 


Using  equations  A-129,  A-130,  A-136,  and  A-l 39,  and  the  orthogonal  itv 


of  the  P ( c os O ) , we  deduce  that 
n 


(x)  - b c'(x)  = md  f'Cv). 
n n n n n 


( A-l 40) 


Dividing  the  right  side  of  equation  A-133  by  the  right  sii.  of 
equation  A-135,  we  obtain 


<l>  (x)  - a 4 (x) 
n n n 

<l>'  (x)  - a C ' (x) 
n n n 


mt*>  (y ) 

n " 

•n(y) 


(A-141) 


which  implies  that 


- mtfi  (v)iji' (x) 
n n n n 

(x)  - mip  (v)t.'(x) 
n n n n 
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which  implies  that 


i p'Mt,  (x)  - ip  (x)c'(x)  = c, 
n n n n 


where  c is  independent  of  x. 


Lemma  A-2.  For  all  complex  z,  we  have 


<P'  (zK  (z)  - <p  (z)c'  (z)  = i 

n n n n 


Proposition  A-2.  For  every  positive  integer  n and  For  x art/ 
y defined  by  equations  A-129  ant/  A-130,  we  have 


n ^^(y)r>n(x)  - m^n(y)C^(x) 


n mif'  (y)  - ^ (y)C'  (x) 
n n n 


'roof  o*  Proposition  A-2.  From  equation  A-133  it  follows  that 


\p  (x)  - me  ip  (y)  = a c (x)  , 
n n n n n 


and  from  equation  A-135  it  follows  that 


i^'(x)  - c ii'(y)  = a c'(x) 
n n n n n 


Dividing  the  right  side  of  equation  A-153  by  the  right  side  of 
equation  A-154,  we  deduce  that 


i^n(x)  - mcn<|)n(y)  4n(x) 
*n(x)  “ cX(y)  " ^n(x)  ’ 


; 
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or,  equivalently, 


(x)c'(x)  " me  (y)c'(x) 
n n n n n 

= c (x)i|>'(x)  - c 4>'(y)c  (x)  . 
n n n n n 

Thus 

^(x)Cn(x>  ~ ^xK^Cx) 

Cn  _ ^(y)?n(x)  - i#n(y)^(x) 

By  Lemma  A-2,  equation  A-157  implies  equation  A-151. 
Equation  A-138  leads  to 


(A-156) 


(A-157) 


*n(x)  - Vn(y)  = Vn 


(x)  , 


and  equation  A-1A0  leads  to 


(A-158) 


'l' ' (x)  - md  ii'  (y)  = b c'(x)  . (A-159) 

n n n n n 

Division  of  the  right  side  of  equation  A-158  by  the  right  side  of 
equation  A-159  results  in  the  expression 


<P  (x)  - d * (y)  A (x) 
n n n _ n 

i)j'(x)  - md  4'  (y)  C'(x) 

n n n n 

Thus 

ip'  (x)c  (x)  - i/i  (x)c' (x) 

, _ _n n n n 

n miJ/'Cy)?  (x)  - i|>  (y)c'(x) 
n n n n 


(A-lbOl 


UV-161) 


A3 


r 


i 


We  obtain  equation  A-152  by  using  equation  A-161  and  Lemma  A-2.  This 
completes  the  proof  of  Proposition  A-2. 

Proof  of  Lemma  A-2.  Observe  that 


j (2)  = [ V — L-.l)N.Z/2)k  ](z/2)v 

v r(k+i)r(k+v+i)JU/i;  • 

K O 


(A-162) 


Thus,  definition  A-lll  implies  that 


ii>  (z)  = (nz/2)*5  7 — il/Al fz/21n('7/'M ^ 

" ' £ r(k+l)r(k-HvM5+l)(  ^ 


k=o 


=v^  i -tv V2> 

V k=o  ^(k+l) r (k+n+3/2) 


(z/2) 


n+1 


(A-163) 


and 


*'(z)  =V7"y  (-l)k(n+k-H)(Z/2)kZn/2n+1 

n V k=0  r(k+l)r (k+n+3/2) 


(A-164  ) 


In  addition,  definitions  A-113  and  A-116  imply  that 


Is  J Jju(z)cos((n+?5)n)-J  , , . 

C (z)  = (ttz/2)4{J  (Z)  - I Stll 

n n+V  ’ 11  sin((n+1i)7r) 


(z) 


] > 

(A-165) 


or 


Vz> 


<7rz/2)  [>Wz) 


l<-»n+V.(lHij)<«>] 


(A-166) 


.4 


1 


Hence,  equation  A-166  implies  that 


Sn(Z) 


-VTI 

k=o 


(-l)k(z/2)k 

r(k+l)T(k+n+3/2) 


(z/2) 


n+1 


+ i 

k=o 


(-11k(z/2)k 

F(k+l)r[k-(n+4)+l] 


(z/2)' 


(A-J67) 


c;(z) 


(-l)k(n+k+l)(Z/2)kzn/2n+1 
* , £ r (k+l)  r (k+n+3/2) 


+ 


i(-DnV^~  I 


k=o 


(-l)k(k-n)2n(z/2)kz  n 1 
r (k+l) r [k- (n+4)+l ] 


(A-168) 


The  nonzero  terms  are 


*n(z)Ci(z)lz-o 


i(-l)nir2n(-n)(^)n+1 
r(i)r(Ss-n>r(i)r(n+3/2)  ’ 


(A-169) 


_ 1 (-!)"* (4) (n+1) 

n (z>  Cn (z) |z=o  r(l)r(n+3/2)r(4-n)r(l) 


(A-170) 


Finally,  we  observe  that 


t (z)i|i'(z)  - ill  (zK’(z)  = A 
n n n n n 

= 1 (-l)nt (n+4) 

r2(l)r(n+3/2)r(^-n) 
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* ■ 


— 


To  complete  the  computation,  we  observe  that 


r(-*s)  = -2  TPs) 


(A— 172) 


since 


r(-VD  = (-4)  r(-4) 


(A-173) 


follows  from  the  known  relation 


T (z+1)  = z r (z ) . 


(A-174) 


Furthermore,  n = 1 and  equation  A-171  implies  that 


i (-1)  (3/2) 

r2(i)r(i+3/2)r(-l3) 


(A-175) 


Application  of  equation  A-174  to  equation  A-175  results  in 


A _ -i(3/2)ir 

i (3/2)  (b)  r (Jg)  r(Ss)  (-2) 


jTT 

r2  (Is)  ‘ (A-176) 


Now  we  see  that  for  Rez  > 0, 


f(z) 


-1  -t_, 
e dt 


(A-177) 
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since 


r(z+l)  = 


tZe  tdt 


z,  -t 
t de 


-tZe_t  + z-1  -t 

L zt  e dt 

1 o ' o 


Observe  that 


= z r(z) 


r(l)  = e_tdt 

' O 


r(4)  = t'^e_tdt, 


Substituting  t = s in  equation  A-180,  we  see  that 


(l/s)e  2sds 


2 e s ds 


\-J81) 


■V  f ". 


since 


implies  that 


or 


(A-182) 


i2= 


TT  / 2 

o 


rdrd9 


= u/4 


(A-183) 


I = Y^/2  . 


(A-184) 


Thus,  equations  A-l/6  and  A-181  through  A-183  imply  that 


At  = i . (A- 185) 

We  claim  that  A = i for  all  n.  Our  claim  is  established  by 
n 

mathematical  induction  on  n.  Suppose  that  n > 1 and  that  A , = i. 

n-1 

Then 


r(4-n) 


T (3/2-n) 

(Vn) 


(A- 186) 


and 


f (n+3/2)  = (n+4)r(n+4) 


(A— 187) 
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imply  that 


A 

n 


i(-l)n-n  (n+4) 

(n+is)  r (n+4)  I'j^yn-) 


Equation  A-188  and  the  original  definition  in  equation  A- 
that 


A - A = U-l)n~1TT(n-l^) 
n n-1  T (n-1+3/2)  f (Ij-(n-l) ) 


This  completes  the  proof  of  Lemma  A~2, 

Previous  results  were  developed  on  the  assumption  o 
form  for  the  expansion  of  a plane  electromagnetic  wave, 
this  expansion.  Each  component  of  a plane  wave  is  a 

the  scalar  Helmholtz  equation 


V2^  + = 0. 


in  spherical  coordinates,  equation  A-190  takes  the  form 


13,23*,,  1 3 

T ?T(r  + -5 ^(sin6  an> 


2 3r 


2 . . 30 

r sin9 


36 


tS  ^ * >'2*  ‘ »• 

r sin  0 3<t> 


Let  us  write 


4)  = R(r)0(0)<D(4>), 


iA-188) 


71  imply 


(A- 18 9) 


f .1  particular 
We  now  verify 
solution  of 


(A-190) 


(t  ) 9i 


• a-192) 
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= l 


a 

m , n , 


pV>ve>y«- 


CA-193) 


and  observe  that  R(r)  = Rm(r),  9(9)  = 0^(0),  and  ( 4> ) = 4>  (0) 

imply 


R(r) 


2 dr(r  R (r))^  + 0(0) 


1 2 . . de 

r sm6 


(sin00'  (0) ) 1 


+ ^fer[~21^*"^)1+k2=  °-  (A~194) 

r sin  0 


Suppose  that  electric  vector  E is  given  by 


E = a^exp(-ikz)exp(iwt) , 


(A-195) 


where  a is  the  unit  base  vector  directed  along  the  x-axis  of  a 
nc 

Cartesian  coordinate  svstem.  The  expression  for  a in  terms  of 

—x 

the  unit  base  vectors  e , e_,  and  e,  in  the  spherical  coordinate 

— r —0  — 9 

system  is 


a = (a  • e )e  + (a  • en.e.  + (a  • e,)e, 

— x —x  — r — r —x  — 0)— 0 — x ^ ^i>  (A-196) 


which,  upon  replacement  of  the  inner  products,  becomes 


a = sinGcosiJ)  e + cos0cos<!>  e_  - sin<t  ei  (A-197) 

— x — r —9  — d>  . 


Substituting  equation  A-197  into  equation  A-195,  with  z replaced  by 
r cos9  , yields 


E = sin9cos(J>exp(-ikrcos0)exp(ia)t)e 


+ cos9cos<t>exp  (-ikrcos0)exp  (iwt)e. 

— V 


- sin<fexp(-ircos9)exp(i<i)t)e^. 


(A-198) 
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Through  use  of  E = M + i N and  equations  A-87  through  A-90 
— —v  — u 

and  A-93  through  A-96,  the  vector  E can  be  expressed  as 


E 


t Og^. 


+ ( 


1 

sin0 


x 3v 

5 3<t>  -0 


3v  , 
30 


. . r 1 rl  3 , 2 3u.  2,  2 , 

+ l {— r-  i—  — (r  — ) + m k ruje 
mk  r 3r  3r  — i 


1 32(ru)  , 1 32(ru) 

— —————  0 ■ 1 - . . - — 0 J , 

mkr  3r30  mkrsinG  3r3$ 


(A-199) 


Our  starting  point  in  finding  expressions  for  u and  v is  to 
write  down  the  spherical  harmonic  expansion  of  exp(-ikrcos0) . 

We  set 


exp (-lkrcosG)  = ) a P (cos0)j  (kr),  (A-200) 

u n n n 

n=o 


where  P (cos0)  and  j (kr)  are  the  Legendre  and  spherical  Bessel 
n n 

functions  respectively.  Now  j^(s)  has  an  infinite  series  expansion 


in(s) 


0n  n c (~l)n‘(n+m)  '■  g2m 

s ^ m! (2n+2m+l ) ! 
m=o 


(A-2G1  ) 


from  which  we  obtain 


ds 


j (s) 
n Jn 


2 "jW 

(2n+l ) 


(A-202) 


Multiplying  both 
upon  integrating 


sides  of 
from  0 


equation  A-200  by  (cos0) sin0  , we  obtain, 
to  n with  respect  to  0,  the  result 
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Vn(kr) 


(cos0) sin6d9  = 


exp  (-ikrcose)Pn (cos9 )sin6de.  (A-203) 


Thus,  the  relation 


r 2 2 

P (cos0)sin8d0  = 

I n 

' O 


2n+l 


(A-204) 


implies  that 


rn 

i exp (-ikrcose)P  (cos0)sin9d0. 
I n 

1 o 


(A-205) 


Letting  s = kr  in  equation  A-205,  we  find  that 


a (s)  (-5— rr 

n n 2n+l 


r 

) = | exp (-iscos0)P  (cos6) sin6d8 . 

I n 

' 0 


(A-206) 


Differentiating  both  sides  of  equation  A-206  with  respect  to  s n times, 
setting  s = 0,  and  using  equation  A-202,  we  obtain 


2%!^  _2__ 
n(2n+l>:  2n+l ' 


1 (-icos9)nP  (cos0) sin0d0 . 
I n 

‘ 0 


(A-207) 


Consequently, 


(-r(I  ^2  ^ 

ant2^lj1_(2^fr)  = (-i)r  ! cos""?  (cos9)sined0. 

1 0 


(A-208) 


By  equation  A-208  and  the  fact  that 


’ o 


~n+l  , ,.2 

cos  0P  (cos0)sin9d6  = — 

n (2n+l) 


(A-209) 


L 
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— m 


we  thus  obtain 


a - (-i) (2n+l). 
n 


(A-210) 


Now  the  functions  u and  v were  previously  expanded  in  terms 
of  products  of  the  associated  Legendre  functions  of  the  first  kind  and 
order  one  and  spherical  Bessel  functions,  P^(cos&) j^(kr) . It  thus 

behooves  us  to  represent  exp (-ikrcoso)  in  terms  of  these  functions. 
Differentiation  of  both  sides  of  equation  A-200  with  respect  to  0 
results  in  the  relation 


r (1) 

ikrsinOexp (-ikrcosS)  = - ) a P (cos0)sin0i  (kr) . 

n n u 

n=o 


(A-211) 


Using  equation  A-211  and  the  fact  that 


Pl(x)  =V 1-x2  P(1'(x), 

n n 


( A— 2 12) 


we  deduce  that 


-lkrsinOexp (-ikrcos0)  = ) a F (cos9)i  (kr), 

n n n 

n=o 


(A-213) 


Noting  that 


r2  71  jn(kr)  + 2r  h jn(kr) 
dr 


+ (k2r2-n(n+l))jn<kr)  = 0 


(A-214) 


and  defining  u by 
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L 


u 


= exp  (iu)t  )cos<j> 


I (-i)n 


n=l 


2n+l  pl 
n(n+l)  n 


(cos0)C^(kr)  , 


(A-215) 


shows  that  the  r component  of  is,  by  equations  A-198  through  A-200, 
and  A-213,  given  by 


E^  = cosi>sin6exp(-ikrcosP)exp(iiijt) 


= exp(iuit)cos<{>  7 
n=l 


(-i) " ( 2n+l ) P 


(cose)  ( 

n 


jn(krl 

-ikt 


), 


(A-216) 


since  by  equations  A-199  and  A-215,  we  deduce  mat  E 


is  given  by 


i r 3 , 2 3u. 

k?[yF(r 


2 2 

k“r"u ] 


r n 

= iexp  (iait ) s<f  (-i) 
n=l 


2n+l 
n (n+1 ) 


n 


(cosA) 


n(n+l)  . 
kr  1 


n 


(kr) 


= exp(iwt)cos<J>  ^ 
n-1 


j (kr) 

(-i)n(2n+l)P  (cosO)  (— 

n -ikr 


(A-217) 


which  is  exactly  the  right  side  of  equation  A-216. 

Now  an  analysis  of  the  coefficient  of  e^  in  equation  A-198  is 
in  order.  An  essential  step  toward  this  goal  is  to  expand  the  function 

F(r,®)  = exn (-ikrcos9)cos0 . (A-218) 

Differentiating  both  sides  of  equation  A-200  with  resnect  to  r.  we 
find  that 


-icos8exp(-ikrcos0)  = £ a P (cos0) j ' (kr) . (A-219) 

n=o  n n 

Since  E = + iN^  and  setting  m = 1 , the  coefficient  of  e 

equation  A-199  is  given  by 
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'■  r . 


(A-220) 


1 3v 
sinS  9<> 


• (L  92u 

k 3r  30 


+ 


1 3u. 

kr  39J 


One  might  be  tempted  to  set 


3v 

3?  ” U 

ar.o  use  equation  A-215  to  deduce  that 


v 


CO 

expt  (iu>t)sin0  £ (-i) 
n=l 


2n+l  pl 
n(n+l)  n 


(cos0) j (kr) . 


Also  one  might  attempt  to  prove  that 


1 

sin6 


•v 

3<i> 


+ i(T 


' r 30 


_1_ 

kr 


) = exp  (iait ) sin0cos<{i 


• exp (-ikrcos6) 


by  substitution  of  series.  However,  there  is  a much  easier 
This  involves  seeking  an  expression  for  u in  the  form 


B (9  , <j>)  exp  (-ikrcosB) 
r 

Then 


~ = B(0,if>)exp(-ikrcos9)  (-  - -kc-°S^)  > 

r 


2 

r“ 


3ii 

3r 


B(0,<j>)exp(-ikrcos9)  (-l-ikrcos6)  , 


J3 

3r 


p 

(r 


3r 


B(0,<f)exp(-ikrcos9)  [ (-ikcosO)  (-1-ikrcosO) 
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(A-221 ) 
(A-222) 

(A-223) 
method . 

(A-224) 

(A-225) 

(A-226) 


- ikcos0 ] . 


(A-227 ) 


Equation.*  A-217,  A-224, 


and  A-227  lead  to 


P _ i ,3  ,2  3u.  .2  2, 

Er  - k7^(r  37)  + k r u] 


- B (0 , <p)exp  (-ikrcosQ)  • {¥~~  [ (-l-ikrcos6)  + 1]} 

+ ikr  [ B (O , <J>)exp  (-ikrcosS) /r  ] 

2 

- ik(l-cos  9)B(0,(!))exp(-ikrcos0) 

2 

= iksin  0B(9,<J>)exp(-ikrcos0)  . (A-228) 


Now  we  desire  to  find  B(e,4>)  so  that  equation  A-228  becomes 


i r 3 / 2 5u  . ,22  , 

T~h- -(r  TT  + k r u] 


kr  3r 


dr 


= exp  (ioot ) sinScos^iexp  (-krcosO)  . 


(A-229) 


The  choice 


B (9 , <f>) 


cosp 

iksin0 


exp  ( iuit ) 


(A-230) 


will  serve  our  purpose.  Thus,  we  must  attempt  to  expand 


, . exp (~ikrcos9) 

u - exp  ( uot ) cost>  ~rr~ 4 

tkrsinS 


(A-231) 


The  method  of  undetermined  coefficients  will  enable  us  to  obtain 
equation  A-215  by  using  equation  A-22°. 

Let 


v = exp(iwt)sin<t> 


exp (- IkrcosQ ) 
ikrsinB 


(A-232) 
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i^et  a(t)  - exp(iait),  and  form  the  derivatives 


TT  * a (t ) exp  (-ikrcosQ)  - ■ C°^  ■ 

ikraine  * 


(A 


= a(t)exp(-ikrcos0)cos$( + 1)  , 

ikrsin^e 


(A 


32  2 
Yr'zs  ~ a(t)exp(-ikrcos9)cos<j>(— <-c°-  j: 

irksin^ 


ikcos6 


+ - c-°-se  ) 

..  2 . 2'  ' 
xkr  sin  0 


(A 


Construct  the  left  member  of  equation  A— 223  by  using  equations 
through  A-235  and  find  that 


1 — + j A 3 u , _1_  3U' 
sin9  dp  Vk  3r30  kr  36^ 


2 

= a(t)exp(-ikrcos6)cos({)( — + - — - s„.° 

ikrsin^e  krsin  6 


+ cos9  + 


COS0 


cose  + 1 ^ 


k2r2sin20  k2r2sin20  kr 


2 

a(t)exp(-ikrcos0)cos<j>( j-  + — c-°-8  ^ 

ikrsin  9 krsin  0 


, isin20 

+ y~  + cos  6) 

krsin  6 


a(t)co30cosil>exp(-ikrcos0)=  E 

0 


(A- 


-233) 

-234) 

-235) 

A-233 


236) 
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I 


This  establishes  the  validity  of  equation  A-223.  Finally,  to  complete 
the  argument,  we  must  show  that 


9v 

3$ 


krsinS  ~a^T<i(ru)  = -exp (iiot ) sinifiexp (-ikrcosQ) , (A-237) 


We  have  from  equation  A-232  that 


= exp(iut)sin$exp(-ikrcos9)  ( — - + 1) 

do  , , 

xkrsin  0 


(A- 2 38) 


Also  by  equation  A-231, 


7"  (ru)  = exp(iwt)sin<))exp(-ikrcos0)  (-— 4) 
ord(p  sinb 


(A-239) 


Thus,  equation  A-239  implies  that 


krsin0  TrTi  = exP(iu)t)sin'»>exP(-ikrcos9)  0 


COS0 

2 

krsin  0 


-) 


(A-240) 


Consequently,  combining  equations  A-238  and  A-240,  we  deduce  that 

9v  i 32 

~ 39  + = -«P(i“08in*exp(-ikrco8e) 


(1 


COS0 


icos9 


-) 


ikrsin  6 krsin  9 


= -exp  (iuit  )sin*exo (-ikrcosQ)  = f • (A-241) 

♦ 

This  completes  the  derivation.  The  following  lemma  has  been  proved: 
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f 


1 


Lerma  A-3.^  If  u and  v are  defined  by  equations  A-231  and  A-232, 
respectively , then  u satisfies  equation  A-2X5,  v satisfies  equation 
A_ 222,  and  u and  v satisfy  equation  A-199,  vector  E is 

defined  by  equation  A-198.  ~ 

In  proving  that  the  u defined  by  equation  A-231  has  the  repre- 
sentation A-215,  we  used  some  properties  of  Legendre  polynomials  which 
are  consequences  of  the  following  lemma: 

Lemma  A-4.  For  every  positive  integer  n, 


2n+1(n:)2 


j^cos  0Pn(cose)sined9  =(2n+1). 


(A-242) 


OKl 


r 2 2 

[P  (cosO)]  sin0d8  = • 

I n 2n+l 

; o 


(A-243) 


Proof  of  Lemma  A-4.  We  use  Rodrigues'  formula 


D , . 1 d,  2 . . n 

P (x)  (x  -1) 

n 0n  , n 

2 n.  dx 

and  observe  that  the  substitution  x = cos9  implies  that 


cosn9P  (cos9)sin9d9  = — - — 
n'  n , 

o 2 n. 


n ,n 

cosn9P  (cos9)sin9d0  =' — — 
„ n -,n  , 
o 2 n: 


= (-D 

„n 


= 

COS0 

(-1 

,n 

n d 

x — 

-1 

dx1 

(-1 

,n 

,d 

( X 

J-l 

, n 
dx 

rl 

(x-l)' 

^-1 

-(x  -1)  dx 


n , . , . n , 


(A-244) 


(A-245) 
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Continuing  the  integration  by  parts  of  the  right 
we  conclude  that 


side  of  equation  A-245, 


it  2n 

cosn0Pn(Cos0)sin0de  = tl) 

0 2 (n+1) (n+2) . . . (n+n)  ^-l'dx 


1 /"a0!  „ 2n 

](x+l)dx 


,2n+l (n ; ) 2 
!n(2n+l).' 


= 2n+1(n.')2 
(2n+l) ] 

To  complete  the  proof  of  the  lemma,  observe  that 


(A-246) 


fPn(cos0)]2sin0d0  = | [P  (x)]2dx 


J_1  n 

1 

f1 

fd^_ 

2 2n ( n ! ) 2 

J_ 

l a n 
1 dx 

(~l)n 

f1 

j2n 

rd  - 

22n(n.-)2 

i-1 

, 2n 
dx 

2 , . n , 


dx 


(x"-l)  ] (x  “1 ) dx.  (A-247 ) 


Since 


j2n 

-1>n  = <2"); 
dx 


(A-248) 


equation  A— 245  implies  that 
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■v  r 


(A-249) 


.2  (-1  ) ' 1 f ' r l ■)  '(  ">  n 

[I*  fensO)]  sinOdO  = ---■ — (x  -t)  dx. 

n 2n,  ,.2 

o 2 ( n I ) - 1 


But  another  integration  bv  parts  shows  that 


fl  0 

(X  -1) ndX  = 

J-1 


(l 

(x-1 )n (x+l)ndx 

J-l 


2n+i(n:r 

(2n+l ) I 


(A-250) 


Substituting  equation  A-250  into  equation  A-249,  we  obtain  equation 
A-243. 


APPENDIX  B 


SAMPLE  PROBLEMS  WITH  COMPUTER  RESULTS 


SABPLE  PR06LEB  1 DECK  SETUP 


CASE  1 (CONTROL  PA8ASETEBS) 

1.000E3  5.986E1  1.0015250  1.00020  1.C00E1  O.OOOEO 


CARDS  2 - (NOC  + 1 ) (DATA  CARDS) 


8. 66C2-01 

5.874E+01 

8.5002*01 

1 • 658E+00 

2.528E+01 

8 . 5005*01 

:. 5995+00 

1 .5795*01 

8.5005*01 

3. 5712*00 

1 . 182E*01 

8.5005*01 

+ . 65SL+00 

8. 9302*00 

8.5005*01 

1 . bSat  + DO 

7.245E+01 

7 . 1572*01 

2.1 79E *00 

8 .651E+0 1 

7.1575+01 

2. 958E+G0 

3.2  31E*0 1 

7.1575*01 

3 . 88  5*CC 

2.8375*01 

7.1575+01 

+ . 7705+00 

' . 365*01 

7.1575*01 

2.598E+CC 

7.  905*0 1 

7 . 6695*01 

2. 95d5+00 

5. 953E+C 1 

7.8695*01 

3 • 571 5 + 00 

8.5565+01 

7.8695*01 

8 ■ 330E  + 00 

3.6075+01 

7.6695*01 

3. 571 E + GO 

8.1955*01 

8. 1875*01 

3. 68 15*06 

6.7015*01 

8. 187E+01 

8, 3305*00 

5.878E+01 

8. 1675*01 

8. 5755+00 

4.5295*01 

8.1875*01 

+.5555+00 

8.3705*01 

8.3665*01 

8. 770  £ + CO 

7.1675*01 

6.3665*01 

1. 6585*00 

7.2855*01 

1.8435*01 

2. 1795*70 

8.6515*01 

1 . 883E+01 

2. 9565*00 

3.2315*01 

1.8435*01 

3. 8815*G0 

2.8315*01 

1.6435*01 

8. 7705*00 

1 .936E+01 

1 .6435*01 

2.1795*00 

7.o745*G1 

8.5002*01 

2. 5985*00 

5.878S+07 

8.5002*0- 

3.279E*C0 

8 .0325*01 

9.500E+01 

8.093E*00 

3.1225*01 

8 • 5005*01 

8.975^*00 

2.5^+E+OI 

8 . 5002*01 

2.998l»C  r 

3.027E+01 

5.9" +5+01 

3.271  ^*  ‘ : 

b. 2775  + 0’ 

5.9045*31 

3 . o + ' Z*  C 

8.4395*01 

5.  085+01 

8.ss3e*;c 

. 9 * 95*01 

5.90+t+OI 

i - . ♦ *, ' 

8.2522*01 

6.  805*01 

8 . 0935+00 

. 505+01 

o..  302+01 

8 . 555*20 

5. O71E+07 

o. 6305+0' 

* . 77  "e ♦;0 

6.3985+ 0‘ 

7. -572*0' 

^’.E'6-*Dr 

7.285E+C1 

7.1575+01 

2. 5 9 a 2 ♦ 0 0 

'‘.6905  + 01 

1.1315+01 

2.951  E+00 

5. 8532*07 

1.1315+01 

7 . 5~1 E+GO 

+.5565+0" 

1.1315+01 

32'03»00 

3. -075*01 

1.1315* 

2 • 956  5*00 

6.0275*01 

3.0965*01 

3.2795*00 

6.2775*01 

3.0965*01 

3.8815*00 

8.9392*01 

3 . 0965*01 

8. 5555*00 

3.979E*01 

3.0962*01 

3.5715*00 

8.1 955*07 

8.5002*01 

3. 88 1 5*00 

6.7015*01 

8.5005*01 

- . 3325*00 

5. +745*01 

8.5002*01 

-. 3305*00 

8.337E*01 

5.8862*01 

8.5555*00 

7.0775*01 

5. 8862+01 

3.57*e*C0 

b. 1955*01 

8.1 30E*00 

3. 881 7*00 

6.7015*0' 

8.1305+00 

-. 3305*00 

5.8785*01 

6. ‘ 305*00 

8.9755*00 

8.5295*01 

6.1305*00 

3.  81E»:0 

8.2525*01 

2. 320E+01 

8.0935*00 

6.6505+01 
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b . 3«05*C0 

8 . 77  0E*O0 

7.  '.675*0' 

f . 3UCB+C0 

8. 7705*00 
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8.975E+00 

7.2855*01 

1 . 643E*0 1 
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The  plot  of  the  computed  absorbed  power  densities  at  internal 
points  of  the  sphere,  Figu-e  3-1,  is  similar  to  that  in  Kritikos 
(2.P-58)  for  a sphere  of  radius  5 cm  and  1000  MHz.  The  elect-ica* 
properties  e and  o are  those  of  biological  tissue  of  high  water 
content . 


Figure  B-l.  Distribution  of  the  power  density  inside 
the  sphere  along  the  z axis. 


t\i  Till 


COPY 


’’POCPA* 

^qjj  O?  PO^’0  I '.'  f.  I * v^»':r'’v’,An-  SPHERE 
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C0M10V  /GJ.TV  V 
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10  50PSXT  fS^O.  ,22  ,?.. 

D0flSGA=:  . D**  ■> .:’’•■ ’ ?c 
SIGHA=SIG*,T^,’:0 

A *= CD 393 ? '?CU?LX  (rPS,  . DL*P — *<*"-.»•'  /’'*  *-  •*  ' ' 

A P G » DO  “ 2C  A ♦ 71  i* r*  * t *.  - 1*  . 

C^-’OC^'.7  -D'"'  ( A ? G ) v (\?.C)  ) 

A1AH0A*2.  < 7924562r'V’  ? 

.\-.T  = 2.7  1 - • :?/*.v  t 

E0= 5C  * . 3 3 3 3 3 ? . . 3 3 2 ? ’* 

X = PI5*D/\L\SD.\ 

Z=AM*X 

C **•  GSV5PATX  SPH^PXCAL  Q5SS2L  v,JVC’*Tr’v~  jvfWT'  •»* *>  wr  (ijjj 

CALL  BESS  (AJ?,Z,N,CQ,S?OPR) 

C *♦*  G^S^PACE  SPH S3ICAL  35SSEL  FDKCTIOVS  J •'*)  » r HV(7) 

DP  (1)  =-DC0S  (X)/X  «- 

DP (2) =DP (1) /X-DSlf (X) /X  > T ' ■ 

DO  " 5 3"  2 r I*  r i" 0 1 ? 2 

AN2  = 2»B-'.  *T10:33 

DP  (S  + 1)  *;.52/X*DP  (F)  -DP  'H-1)  IE  T :t 

IF  (DAOS  (DP  (Li))  .GE.  STOP?)  GO  TO  2?  "1*0035 

15  COtilIVCS  BISC'?6 

20  N P 1 = B ♦ 1 51*2037 

n=b  ri-ooje 

P (BP*  ) =>C  . DO  HI50039 

? (hpi-1)  =-i . dc/  (x*x*d?  (»?i) ) o: 

NB2  = NP1-2  *1*704* 

DO  25  B*1,KB2  *1*0042 

I*VP1-B  «iEC',a? 

AV2=2*I-1  SIEOOut 

25  P(I-1)=AK2/I*P(I)-P(I*1)  KIE0C4S 

FQ=DSIS  (X)  / (X«p  (1)  ) *157745 

DO  3fl  B*1, HP1  *.150247 

30  P (B)  =B2*P  (fl)  «'IEC0#e 

N*1=V-1  *150040 

K*2=V-2  »I5C''5C 

C0*1. DC/CO  * 1500  51 

RC  ~ 1 • DC/PQ  *150752 

0 **»  P? 1ST  001  15,  PA5IC  I*r"T  DATA,  AVD  EFFO®  AHALTSIS  p\’(  *1*005? 

PHI ** T 35,EESQ,ALlBDA,50,5IG(n ,EPS#D,TXHE# AR,IB1,X,P0,CC  BIE005. 

35  POSBRT  ( *1 DEPOSXTIOB  OP  OO-DH  IVOIOE  A BOBOGEKEO!  SPHERE  XBSBBSF  IE  55 

ID  IV  A v BtECTFOBAGlETIC  FI EID •/ > - *"000 - v «■,*?.  m V\"-' 

*LSRSin  5,’  CB  FIELD  SIS?  »/B,/,0"’V  OC7BI! 

Ilf ITT  »•,.  * /•  IBLATIf  DIELECTRIC  'ORS1AVT  •• ,F7. 2, • BIE0058 

1 3i  • CB  ~l*3  * • , F7. 2 , SEC/'OREFR-  *1*0059 

1IHPSX  ••>1P2E13. 5/*  T Z HIGHEST  R B»IH  * CTIOi  CO ■ LSRIE0960 


tJ  o 


1 OF  ORDER  -',13,'  SIZE  PABARFTEB  - • ,E1 3.  5/' C 8 ATIO  OP  THE 

* CALCOLATIOa  OP  THE  ZEBO  ORDER  SHE8ICAL  BESSEL  POBCTIOR  OP  A6 
IT  TO  THE  VALUE  OF  SIB(I)/X  = • , El  5. 5/'  OBATIO  OF  THE  PIRST 

10LATIGN  OF  THE  ZEFO  ORDER  SSEBICAL  BESSEL  POBCTIOR  OP  ABGUREB 
1*  TO  THE  VALUE  OF  SIIi(Z)/Z  **,2F15.5//) 

C **»  GEHEBaTS  COErPICIEHTS  C*  AND  DR 

DO  nn*1,NR2 
NP1=NN*1 
bpf.x=i*p  ;«pi) 

FPSI=DCRPLI  (RPHI , -E*DP  (NP1) ) 

CPHT=Z»AJB (SPl) 

A N = N N 

F A C T = 1 .DO/ (2.  CO*  AN*1  .DC) 

P PH  I P =P  'PI ) *X»  (FACT*  (AN*P  (BP1-1 ) - (AS*1.  DO)  *P  (NP1  *1)  ) ) 

RaIP  = -.'P(S?1)-X*  :-ACT»  (AH*DP(RP1-1)-<A»*1.D0)*DP(BP1*1))) 

F?SIP*DCftPLX  (RPHIP.RKIP) 

cph:p=ajr  (npi , *z*  ifict*  , »n*ajb  (npi-1) - <ar*i.do) *ajb  (*pi*1) ) ) 

Q-  ( B PHI  P*  BPS  I “BP  HI*  EPS  IP) 

CN  < V N ) =0/  (CPHIP*S-r'SI-Aft*CPBI»RPSIP) 

43  DN (NN) = Q/  (AS»C?HIP*aPSI“CPHI*BPSIP) 

IF  (NOC.  EQ.  0)  GO  TO  65 
DO  60  0=1, HOC 
READ  10,  B , THETAD , PHID 
D1=2. DO*  P 
FHI=PHID/RAD 

IF  (THETAD. EQ.OLDTH)  GO  TO  45 

THE7A=THETAD/BAD 

SIRTH=DSIB (THETA! 

CALL  PL  (THETA  ,N,?,DP) 

45  IF  (D1 . EQ. OLDD1)  GO  TO  50 
RKP  = ALAHDA/  (PIE»D1) 

Z=  A ft"  PIE*D1/AL»SDA 

CALL  BESS  (AJR,Z,RC,Q,  STOPS) 

NC=  NC-2 

IF  (MC.GT. NK2) BC=NR2 

C ABSORBED  PONER  DENSITY  AT  GIVES  IHTKRMAL  POINT  OF  SPREB 

5C  CALL  EVEC(PD) 

PD=.r5D^*SIGHA*PO 

C •«»  PRIST  COT  INTERIOR  POINT  PABTICOLABS 

P ? I P T 55, B ,Tr.2TAD,PnID,?D 

55  FOR  ft AT  (‘  INTERNA.  POINT:  RADIOS  *«,P8.3,*  Cft  THETA  = • , F7 

-DEG  PHI  = 1 , ?7. 2,  ' DEG  ABSORBED  PONEB  DENSITY  *«,P12.8, 

1TS/V  *5*) 

I CONTI  NOE 

65  IF  (I OPT.  EQ. C ! GO  TO  5 

TOTAL  ABSORBED  PONES  AND  AVEBAGE  ABSORBED  PGR  ER  DENSITY 
TOTPOK  = .'!5D-f.«SIGRi*GADSS3  (83) 

PAVG=?OTPOB* 1 ■ D6/ ((4.D0/3. DO) *PIE*  (D/2. DO ) **3) 

•***  PB  1ST  OUT  ABSORBED  TOTAL  POBEB  AMD  AVEBAGE  ABSORBED 

POKES  DENSITY 
PRIST  75,TOTPON,PAVG 

75  FOR  ft  A T CO  • ,91, 'TOTAL  ABSORBED  POBEB*  • , 1 PEI  3.  5 , • BATTS* /' 0 • ,9 
1 ER  AGE  ABSORBED  POBEB  DE NS ITY • • , E 1 3 . 5 , • BATTS/H** 3 • ) 

GO  TO  5 
80  STOP 
END 

SUBROUTINE  EVEC(PD) 
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CALL  T7R«  (>-C  , V) 

X2=T71*OrPJv * 

CALL  T2?B  (HCK#Tr) 

URE=Dr^*T2 

T3«TS*  (-2.  V * BES2H1/I-  (1 . D?-”AC2/ZSQ' 

CALL  TERM  (TCP  ,T3) 

0PB*«u3?S+TJ 

I?  (CDA33  (T'>  .GT.CDA8S  (UE)*S?P)  GO  TO  5 
IF  (CDA3S  "T  .(>?.C*'ASS  (?!,S**E?S)  CO  TO  5 
IP  (CDASS  (??)  . T.E.CDXBS  iaERB)*BPS)  IE=1 
5 IP  (IT.  10. 1)  GO  TO  45 

IF  < (THr  A.  GS.  1.  OD-6)  . APO.  (THETA.  IT.  2.  141592D01  ) GO  TO  20 
"‘T*DN  (P<  "'Cl/2.  DC 

if  ?.  * u * ?92d0)  t?=  (-1  • 7 “ * *•  vvpi*  t 

GO  TO  25 

2C  TTSDV  ' •"  ) * 9KC3*  ? (»*v>  /c x v~'j 
25  T1  =?T*  4*  ’ ' * •*■  .) 

CALL  ‘ ?«  (:  CK,?1) 

0?  = u?f-, 

TT«CM  O ,)*  "AC3*!)P  (PN) 

?2  = ?T‘*  *.J 9 fvV*1) 

CALL  ?H?.,r  . icr ,72) 

usr*tm*r2 

T3=TT*OI  VT 
CALL  TEPR f ' CK #T3) 
tJ?PT*l’S 

IF  (C  . 3 f (7 * ' . G? , Cr n ^ f ?S)  CO  TO  35 
IF  (C0A3S  (T2)  .CT.  CD1  ' . (U  " -*7PS)  rO  *0  '5 

IF  (CLAPS  (T2 j .-V.C?A5S(UH  ?) *2P3)  IT*1 
35  I?  GO  TO 

TpsDK  ' . *PAC3# D» *v  > 
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* IF 01 77 
RIF0178 
* 'S0179 
“IE0190 
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<.Kk.'  h : a scfc.T"*) 

u?*i . ♦'" 

If  iiTH&Tk. GE. 1.CD-6) .AND. (THETA. LT« 3, 141592D0) ) GO  TO  6C 
7P»C K (NN) * PaCI/2. DO 

If  *7  BE?  A.  GE.  3.  1 41 592  DC)  TP  = (-  1 . DO  ) * • MN  P 1*  TP 

GO  TO  Sl 

oC  7?*C*‘  NN  *Pf.C34P(NS)/SINTB 
65  ?2*  TP*  . i *<  N*1 ' 

C A L ».  TERR  .i  C K , 7 * j 

U K P - 

T 3 = TP*HRJN' 

C.'LL  T :■:.<«  (NCK,?3) 

RP^'JHP  ? .1 

•*  r"  . .T,cOAt>:  (up}  «e?s;  go  to  6a 

7 ^ • : ■:>  . G7.CDABS  (0RP)  *Er3>  GO  TO  6 6 

; c 3)  . , E.  -0/B5  ( JRRP)*E?S)  IP«1 
c "♦I*',  iw.  3>  oO  TO  71 

<n^jr=  ac..* 

7 . F ( < C K . GT.  4 N C A * 1 

7*  /C 0 = £,'  " . PHI) 
iz  = :. 

2f-r:.*7.  ■ • . . 0 *Aft40F3*AHZ*  (URRE^OB)  ) 

*.■  -EO^DCAPH  1-7  . i AG  (ER)  , DBEAL  (ER)  ) 

Sl»9*C.:X  * FCOw  (BKR*CRT*AM40P*7) 

••TH'J?.  =-C*  A*C2X4FC0*U74DCrtPLX  (-DIAAG  (SUB ) ,DREAL  (SDK)  ) ) 

ESI *0 SIN  (Phi) 

30ft*-C2X*?SI4  (RXJ 4nR?*A.‘.*ORRP) 

EPHI*BO*  -AH*CEX*  I*UP«*DCHPLX  (-DIBAG  H)  0 X JH)  ) 
r 0 = OF  SAL  (2R* DCGNJG  (AR) ) ♦DBEAL ( ETH2T A*OCONuG  (ETHE7A; } ♦ORE A L (2  PH  I 
3 ON JG  (EPHI) ) 

RETURN 

END 

SUBROUTINE  TERN  i A C X , ? ) 

COR PG7ES  {- J) 44  N*  ( N-TH  TEBfl  IN  SERIES) 

IMPLICIT  REAL4 6 (A-H,0-Z) 

COMPLXXMt  T 
GO  70  (5,: 0,15, 20)  ,MCK 
5 T -DC*  PI  I (DIRAG  (?)  ,-DBEAL  (T)  ) 

GO  70  2^ 

If.  T=-7 

GO  7 •"  2* 

1-5  ^-HAaG  (?)  , DR2AL  (T)  ) 


• * „0j'.  . ..  s v -o  . . , y , stop? 

C . .0  SPhlBICJ-L  bESSEi,  FUNCTIONS  ok  (BX)  AND  NN  (RXj 

*'  -/•**  (a-H  , 0“Z) 

A*.*  / CO i 

■'■  P _ — i . t. , , K f YO  , Yl  f Q 

-~c;  ,vS(L*/z 

-CDS I h 'Li 

s (Yr,-Q,  /Z 
00  5 .1**  , uc 
7=  (2*!  ft-?  i /O’  Y1-X0 

IF  (CDAuS  (Y  ).G2.  STOPS)  GO  TO  10 

Y“*  V. 

6 Y 1 = Y 

*C  IF  ( A . GT.  3)  GO  TO  jO 

C *4*  PRINT  OUT  EBROB  RESSAGE 


o o 


{ 
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PRINT  25, Z 

25  FORMAT  co*»»*  PROCESS  CAR  NO*  PROCEED  SINCE  RB2 
1 * 5.  71 

STOP 

30  AJ  (Ml  *DCMPLI  (O.DC  ,P.  D91 
A J (N-’'*-'.D0/(Z»  Z*T) 

NM2=N-2 


m:e02ui 

0 POP  Z «• ,’P2D~IS0242 
MIB02«' 
MIE0244 
MTX0245 
NIEC246 
»TZC247 


DO  35  1*1, NM2  FIE0248 

L = K-I  HIED249 

35  AJ  (L-1)  = (2*5-1)/Z*lJ  (1) -AJ(L»1)  *1X0250 

C*0/  (Z*A.T  ('  • ) « ISO 2 SI 

NM1=M-1  MIS0262 

po  4C  *(-’,<1x1  ri?3:53 

aj  m =o*aj  <m  »:»C254 

F fCT  BS  (AJ  (8)  ) , l*.  s . P-25)  7ETUER  » " ®C  255 

«:  ccn~i '■■•»  *x?0256 

PETUPW  s 1X0257 

"PO  «ir0258 

?U5POCTTRE  fTHETA  ,M,P,D55)  '•  T2025P 


ppv  7BATrS  / S'Tr*  " ,jw;*wpr'P  ’ v r «ur«Tp 

PIPS"  DEPIVA1 . »•••  i 
IMPLICIT  "AL' 8 (A-  ,0-Z) 

^0  5!  T ' ' C' 

‘N\  - DSTN  (*f  T~px) 

C».’=DCOS  (Tl"^7.A) 

P ’ , »F*2 

? .2)  =3.  90»?NO*CP>.T 
op  (i)  =csj 
OC  50  M=2,H 
A = M 

MP’=M*1 


P f"*1)  = (2.  DC*  A*1 . DO) /A«CNJ»P  (»' - (A*->.  DO) /A*P  (5-5)  * ‘E0272 

IF  ( (THETA. GE. 1.CD-6) . AND.  (THETA. LT. ?. 1 ;’$92D0) ) GO  TO  5 M *0273 

P?  (M)  =H*MP1/2  MIE0274 

IF  (TBETA.«E.  3.  141592P0)DP(M)  = (-1 . D0>  **  M*  DP  (M)  MIE0275 

GO  TO  10  MIE02"6 

5 DP  (M)  * (A»CNJ'*P(M)  - (A*1  . DO)  *P  (H-1)  ) /SNJ  MIE0277 

1"  CONTINUE  MIZ0278 

BETP5N  MIE0279 


•IED2! 


* X ' 0 X * c 

* ~ r>n  r 


' '•0268 
*150269 
” I“0270 
* ;"0271 


END  MIZ028C 

FUNCTION  GA0SS3(H3)  MIE0281 

PFO  DUCT  SOLS:  M-POIRT  GAUSS-LEGENDPE  QUADPATUBE  FORMULAS  MIED282 

IMPLICIT  PEAL*S  (A-H,0-7.)  MIE0283 

COMMON  /GAUSS/0 , Ml , M2 ,R  MXE0284 

COM*  ON  CR  (1 r.  P)  ,DN  (IOC'  ,AUS  (5f0)  , A»  ,AMR,CFX,Z,  DP(1?0)  ,P  (If  15  , ALAHPASIE028S 
1 ,D1 , SO, PHI, PIE, BKB, STOP  , IHTH.THE* A,»C,»M2  RISC 

COMPLEX* 1 6 CN ,DN, AJP , AM, AMK.CEX.Z, ) MIS0287 

TIMERS  ION  RPOINT  (91  ,f.Ev  1 0)  , T ( 3 7 ) , JT  ( 3 3)  , A RG2  (2)  , A PG  3 (2 ) MXE0288 

DAT!  NPOIBT/2,3,4,5,6,8,13,12,14/  HIB0J 

DATA  KEY/1 ,2,4,6,9,12,16,21,27,34/  •X’’029C 


D!~A  T/0 . 577  3r0;.'f. 8962DC  , 0.  DCOr  X I "oDOOOf.  DC  , 0. 77oS4«,f,s<>2  c *#o  r>M  (p  . . •<.» 

1 33998 10 4358486 DC' ,C • 861  * 36311594G50  .C  0000  OPCOC  , ,5  8469310M 

1X056  DO, O.90617984593866DO, 0.23861  1 6083201  0.61  ’ 1 662  DO , C' 

1,9  246951420  3 1 500  , 0. 1834396429956!  )0,0. 5255  290991633D0,  . 796666978 IE0294 
179  3D0 ,0. 96P 28985649  5 400 ,0.'» J87433898163O0 ,0.  433395394124250P, BIE03  95 

if  .6794095682990200, 0.865063366689 9800, 0.9739065285171 7|>r>,0. 1252339BIE0296 
If 851 14700 ,0.367831498  }1  00,0.  ’ 954286621  , "<  19439001)3  0297 

,C.  9(  41172S677C970C,P  . 15- 06J42a*>'2"?  , . If  • 1 Xu  )4  •'’O'*  34  DO,  D.  '91128180296 

1 36  3 92’89DA,C . 5152486363581500,0.!  I I 9 481 16900,0. 8272013150*  nbOHIS0299 
10, C. 9284348836635700, 0.9862838086  68100/  MIE0300 


b 





Av 


DATA  «T/1  .CCi'CCP  CCOOOOODC  ,0.88°- 86R888888SI?  , 0.  55*555 555  555b 5 Dv , * .1I203C  1 
1 . 652’ 451 5486255D0 ,C. 3478548451  3"’4  5D0  ,0. 5686b8886888R9DG,C  . 47 06 2867b 12 030 2 
1C  69  4 3 706,0. 2 38  926  8850  56 1900, 0.46  7 91  39  36  57  26900,0.  360761  5730481 4 DO,  «.  ISO  303 
1C.  171  3246923791700  ,0.  36268378337  8 36D0, 0 . 31  370  6645  87789DC  , 0 . 222  3610 8IE0  306 
134  45  3 3700,0.  10 1 228536 290 38D0, 0.2 9 5524224716  7500,0. 26926671 93 1900  DOST  20 30  5 
1 , 21 906636 251 59 8D0, 0.149 45 134915058D0.C.  066671 34430 86 920,0.2491475120 30 6 
1C4581  340  DO,' . 2334925365383500,0. 2031674267230700,0. 16007e32654335DHIE0 307 
10,0. 1067393259953200,0.0471753363865120,0. 21 526 3853 4631 6D0,C. 205198120308 
1 8 46 372130  DO, 7.  1 8 55383 9747794  DO, C . 157203 16 71 58 1900  ,0. 12151 857G6879CBIE0 309 


1DC, 0.0 80 15808715976 DO, 0.0351194603317500/ 

RIE031G 

20  5 1 = 1 , 9 

8IE031 1 

IF  (R3.E0.SPPI*T(I))  go  TO  20 

RIE0312 

5 

CONTINUE 

8IE031 3 

1? 

PRINT  15,  .11  ,82 .83 

8IE0314 

15 

FOB  8 A 7 ■•-3FF0R  IN  INTEGRATION  CONTROLS.  81  =,,I6,1 

82  =• , 16 , 1 

BIEG315 

183  = ■ ,15) 

8IE0316 

GAOS  S 3=0 . DO 

RIE0317 

RETURN 

8IE0318 

20 

JP3  "KEY  (I) 

BIE0319 

JL3  = K£Y  (I  ♦ 1 } -1 

8IE0320 

DO  25  1-1,9 

BIE0321 

IF  , 8 2.  EQ.  NPOINT  (I)  ) GO  TO  30 

RIE0322 

25 

CONTINUE 

8 ISO  323 

GO  TO  10 

8120324 

30 

JF2  "NET  (I) 

8120325 

JL2  = KEY(I*1)-1 

H1EC326 

DO  35  1=1,9 

BIE0327 

IP  (R1.EQ.  NPOIVT(I)  ) GO  TO  40 

8120328 

35 

CONTINUE 

SIF0329 

GO  TO  10 

8120330 

40 

JF1  =KSY(I) 

BIE0331 

JL1  =RBT(I*1)-1 

8120332 

INTEGRATE  0V2S  TKETA 

8120333 

P02=?I2/2. DC 

8120334 

S = D/4.D0 

BIE0335 

SUK  3=C . DO 

8120  336 

DC  85  J3=JF3,JL3 

8120337 

I?  (T  (J3)  . NE.O.DO)  GO  TO  45 

HIE0338 

A3G3 (1, = ?D2 

SIE0  339 

13  = 1 

8IE0 340 

GO  TO  V 

8120341 

45 

ARG3  <1)  = ?D2*PD2*Y  (62 

8 1 EC  342 

ABG3  (2)  = ?02-PD2*Y  (03) 

8120343 

13  = 2 

8120344 

53 

DO  85  L=1, 13 

8I2C  345 

THETA";..  G3  (l 

8120  346 

SINT  H =2SIN  (T’.-.ETA) 

P.  120  347 

C A i.L  ?L'7HETA,N,P,D?) 

8120348 

INTEGRATE  GVE8  RADIUS 

BIE03U9 

3 u 82=' . o' 

SIE0350 

DO  80  JI=JF2,JL2 

SI203S1 

IF  (Y  ;J2:  . NE.C.DO)  GO  TO  55 

8120352 

A?  G 2 (1)  = 8 

8IE0353 

12=  ‘ 

8120 

GO  TO  60. 

8IEC355 

55 

APG2  (1 ) =R*  Y (J2)*R 

P.I2C356 

ASG2  (2)  = R-T  (J2)  *B 

SIE0357 

12  = 2 

8I203''3 

60 

DO  80  1=1,12 

BIID3S9 

01  = 2.  DO*  48G  2 (I, 

8 120  360 

77 


